Botany 2020

Setup

Install

Prior to starting this workshop, you should have had the following packages installed.

## Make a list of packages
list_of_packages <- c("tidyverse", "ridigbio", "spocc", 
                     "scrubr", "RCurl", "rjson", 
                     "sp", "raster", "maptools", 
                     "maps",  "rgdal", "mapproj", 
                     "ggplot2", "ENMeval", "dismo", 
                     "RStoolbox", "hypervolume", "gtools",
                     "caret", "phytools", "factoextra",
                     "FactoMiner", "alphanull", "picante", 
                     "devtools")

## If you do not have these packahes installed, install the package
new.packages <- list_of_packages[!(list_of_packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

## Make sure all packages load
lapply(list_of_packages, require, character.only = TRUE)

# Install packages directly from github
library(devtools)

## ggbiplot
install_github("vqv/ggbiplot")
### If you have R v.3.6
#install_github("vqv/ggbiplot", ref = "experimental")
### If you are having issues installing ggbiplot, you may have to uninstall the package 'backports' and reset R. 

## ENMTools
library(devtools)
install_github("danlwarren/ENMTools")
library(ENMTools)

# If ENMTools does not install 
# library(devtools)
# install_github(repo = "ecospat/ecospat/ecospat", ref="3.1")
# library(ecospat)
# install_github("danlwarren/ENMTools")
## If you are still having issues, you may need to reset R and try again. 

## ENMeval
library(devtools)
install_github("bobmuscarella/ENMeval@master", force = TRUE)
## rSDM
devtools::install_github("Pakillo/rSDM", force = TRUE)
## spThin
library(devtools)
install_github('johnbaums/rmaxent',"spThin")

What are these packages?

tidyverse is a collection of R packages that are used for “everyday data analysis” including ggplot2, dplyr, tidyr, readr, purr, tibble, stringr, forcats. We use dplyr and ggplot2 throughout the many R scripts.

ridigbio is an R package that query the iDigBio API.

spocc query data from GBIF, USGS’ Biodiversity Information Serving Our Nation (BISON), iNaturalist, Berkeley Ecoinformatics Engine, eBird, iDigBio, VertNet, Ocean Biogeographic Information System (OBIS), and Atlas of Living Australia (ALA).

scrubr can be used to clean the data.

RCurl and rjson are used to query an application programming interface (API), or an online database.

Geospatial packages include sp, raster, maptools, maps, rgdal, RStoolbox, and mapproj.

ggplot2 is used to plot and visualize data. ggbiplot is an add on to ggplot2.

ENMTools, ENMeval, and dismo are specific to ecological niche modeling.

hypervolume and caret have functions related to statistics.

phytools and picante are related to phylogenetics.

factoextra and FactoMiner are used for statistics related to a principal components analysis.




Troubleshooting

This exercise required the package ENMTools. This package requires Java (and the Oracle Java JDK) and that the maxent.jar file is in the folder that contains the dismo package (this package is a dependency of ENMTools). To find out where to put the jar file, run the command

system.file("java", package="dismo")

Once you have your java and maxent file set, you are ready to go! If you have more issues, check out: https://github.com/MTFA/CohortEx/wiki/Run-rJava-with-RStudio-under-OSX-10.10,-10.11-(El-Capitan)-or-10.12-(Sierra)

To download the jar file to the right directory —

library(devtools)
install_github("johnbaums/rmaxent")
library(rmaxent)
get_maxent(version = "latest", quiet = FALSE)

Then reload rgdal (it is masked)

If you are having error with rJava

## Remove rJava
remove.packages(rJava)
## Update your R Java configuration 
system("sudo R CMD javareconf")
## Restart R and check to see if the path matches
Sys.getenv("DYLD_FALLBACK_LIBRARY_PATH")
## Reinstall rJava
install.packages("rJava")
library(rJava)

If you are still having issues with rJava - check to see if you are using Java 12 (this should be included in the error message)

system("java -version")

The inital relase of java 12 does not work - install an old version instead.
Install an old version instead Find on the Oracle website or use homebrew.

system("brew cask install homebrew/cask-versions/java11")

Restart R and redo the previous steps. Required: unistall the newer versions of java prior to repeating and restart R after.

## Uninstall Java 12, then repeat the following
## Remove rJava
remove.packages(rJava)
## Update your R Java configuration 
system("sudo R CMD javareconf")
## Restart R and check to see if the path matches
Sys.getenv("DYLD_FALLBACK_LIBRARY_PATH")
## Reinstall rJava
install.packages("rJava")
library(rJava)

If you have additional issues, let us know!

Download Occurrence Data

Original script by AE Melton.
Modified and created by ML Gaynor.

Load Packages

library(ridigbio)
library(spocc)
library(scrubr)

Data download from iDigBio

First, we are searching for the species Liatris cylindracea. Next, download occurrence records for the family Liatris belongs too, Asteraceae.

iDigBio_LC <- idig_search_records(rq=list(scientificname="Liatris cylindracea"))
iDigBio_LC_family <- idig_search_records(rq=list(family="Asteraceae"), limit=1000)

What if you want to read in all the points for a family within an extent?
Hint: Use the iDigBio portal to determine the bounding box for your region of interest.

First, set up your search input. Check out the ridigbio manual to understand the input parameters.

rq_input <- list("scientificname"=list("type"="exists"),
           "family"="asteraceae", 
            geopoint=list(
                    type="geo_bounding_box",
                    top_left=list(lon=-87.86, lat=30.56),
                    bottom_right=list(lon=-79.21, lat= 24.78)
                    )
           )
iDigBio_LC_family_florida <- idig_search_records(rq_input, limit=1000)
uuid occurrenceid catalognumber family genus scientificname country stateprovince geopoint.lon geopoint.lat datecollected collector recordset
006b407f-aafc-49e5-82c9-31ff8aff3fe3 b70c7beb-8a13-4779-8a5d-9e952edc75e9 brit287702 asteraceae liatris liatris cylindracea united states NA NA NA NA NA 72d4c3c7-3413-4588-a803-e1a63e0d7c6c
00cb52c6-3606-47ec-94cd-72222ac77f2c 46d67299-9122-4bb8-80af-1ac000ff76d4 v0033039wis asteraceae liatris liatris cylindracea united states wisconsin NA NA 1958-07-20T00:00:00+00:00 hartley, thomas 9e5aede6-bee5-4a3d-a255-513771b20035
01f04596-0721-4556-8831-86a39cbb2c58 373bea21-5138-4730-8915-2347f84e5256 ind-0137239 asteraceae liatris liatris cylindracea united states indiana NA NA 1915-08-22T00:00:00+00:00 chas. c. deam 0dab1fc7-ca99-456b-9985-76edbac003e0
021b6c1d-9b22-45e4-9cea-a26821e65910 urn:uuid:788e76f1-4e00-4ead-b81a-b30179edff26 211373 asteraceae liatris liatris cylindracea united states minnesota -92.83619 45.42608 NA taylor, b. e2def7e2-1455-4856-9823-6d3738417d24
0237adc6-dd85-4f87-9e5c-422996272a9f a48f306d-d338-471d-a2ee-82100933d41b b99374 asteraceae liatris liatris cylindracea united states wisconsin NA NA 1937-09-08T00:00:00+00:00
  1. pohl
c2767dde-1315-4d78-abf9-8e098dd588ab
039f4bde-2ab1-4561-a9af-1c850a42aab9 urn:uuid:8c635507-6f3d-4936-b642-221666aafb26 211375 asteraceae liatris liatris cylindracea united states minnesota -91.31853 43.62916 1919-08-01T00:00:00+00:00 rosendahl, c. e2def7e2-1455-4856-9823-6d3738417d24

Data download using spocc

Using the package spocc we download data from iDigBio and GBIF. We specified that coordinates are needed. Notice that the occurrence records for iDigBio and GBIF are stored sperately.

(spocc_LC <- occ(query = "Liatris cylindracea", from = c('gbif','idigbio'), has_coords = T))
## Searched: gbif, idigbio
## Occurrences - Found: 579, Returned: 579
## Search type: Scientific
##   gbif: Liatris cylindracea (480)
##   idigbio: Liatris cylindracea (99)

See the spocc manual to set the extent for a search.

Processing spocc download

Synoymns

Sometimes the ‘name’ in the database does not match the ‘scientific name’- this can be due to synoymns or variety names.

##      [,1]                         [,2]                        
## [1,] "Liatris cylindracea Michx." "Liatris cylindracea Michx."

Since this could be an issue when processing the data, we are going to fix the name using the package scurbr.

spocc_LC_namefixed <- fixnames(spocc_LC, how = "query")
## Warning: fixnames will be removed in the next version; see scrubr::fix_names
Data format

Since the occurrence records for iDigBio and GBIF are stored sperately, next we have to convert the occurrence download into a single data.frame (df) where all query are combined.

spocc_LC_df <- occ2df(spocc_LC_namefixed)
name longitude latitude prov date key
Liatris cylindracea -90.62685 38.26096 gbif 2020-03-18 2596076340
Liatris cylindracea -87.16622 33.18885 gbif 2020-06-14 2802751979
Liatris cylindracea -89.49381 42.69785 gbif 2019-07-03 2283149029
Liatris cylindracea -87.39522 41.62323 gbif 2019-07-13 2294490090
Liatris cylindracea -89.77279 43.13965 gbif 2019-07-19 2294629302
Liatris cylindracea -89.39378 43.18260 gbif 2019-07-23 2294645919
Liatris cylindracea -79.42030 43.75006 gbif 2019-07-23 2294650703
Liatris cylindracea -87.89371 42.45455 gbif 2019-07-25 2294679719
Liatris cylindracea -83.09190 45.98512 gbif 2019-07-24 2311314905
Liatris cylindracea -88.20168 42.06006 gbif 2019-07-26 2311338809
Liatris cylindracea -84.49025 44.62769 gbif 2019-07-25 2311387297
Liatris cylindracea -87.82655 42.53831 gbif 2019-07-25 2311413181
Liatris cylindracea -87.80137 42.45064 gbif 2019-07-25 2350326766
Liatris cylindracea -88.31463 42.19416 gbif 2019-07-07 2397483624
Liatris cylindracea -81.93469 46.09624 gbif 2019-07-14 2447982139
Liatris cylindracea -80.04160 43.44048 gbif 2019-07-31 2465114768
Liatris cylindracea -81.81443 43.21843 gbif 2019-07-29 2540734101
Liatris cylindracea -82.94771 45.92069 gbif 2019-08-01 2331855805
Liatris cylindracea -87.80560 42.40589 gbif 2019-08-02 2331867098
Liatris cylindracea -80.20854 42.85761 gbif 2019-08-03 2331914267
Liatris cylindracea -82.77755 45.69730 gbif 2019-08-04 2331944192
Liatris cylindracea -86.03655 42.60887 gbif 2019-08-04 2350326417
Liatris cylindracea -81.82243 42.39630 gbif 2019-08-05 2350326885
Liatris cylindracea -87.39025 41.61772 gbif 2019-08-03 2350327056
Liatris cylindracea -87.39561 41.62389 gbif 2019-08-03 2350327166
Liatris cylindracea -87.54515 41.75903 gbif 2019-08-06 2350328405
Liatris cylindracea -82.69163 45.61977 gbif 2019-08-08 2350331158
Liatris cylindracea -89.43021 43.04116 gbif 2019-08-09 2350338061
Liatris cylindracea -87.27676 41.61042 gbif 2019-08-10 2350353575
Liatris cylindracea -81.94163 43.21775 gbif 2019-08-10 2350354522
Liatris cylindracea -81.97216 43.27360 gbif 2019-08-10 2350354713
Liatris cylindracea -81.92199 43.37881 gbif 2019-08-10 2350355970
Liatris cylindracea -86.29775 35.78874 gbif 2019-08-08 2350356515
Liatris cylindracea -81.87154 43.21321 gbif 2019-08-09 2350359030
Liatris cylindracea -88.25393 42.01357 gbif 2019-08-10 2350361500
Liatris cylindracea -81.90357 43.29497 gbif 2019-08-10 2350373293
Liatris cylindracea -88.20179 42.05989 gbif 2019-08-11 2350377389
Liatris cylindracea -88.49212 41.83114 gbif 2019-08-10 2350384931
Liatris cylindracea -87.55979 41.67936 gbif 2019-08-12 2350389298
Liatris cylindracea -87.55978 41.67948 gbif 2019-08-12 2350391421
Liatris cylindracea -81.83265 42.27649 gbif 2019-08-12 2350404366
Liatris cylindracea -81.97914 42.29385 gbif 2019-08-14 2350423491
Liatris cylindracea -81.85499 42.22794 gbif 2019-08-14 2350437522
Liatris cylindracea -82.68709 45.66322 gbif 2019-08-14 2350439840
Liatris cylindracea -87.90551 42.46278 gbif 2019-08-03 2350444378
Liatris cylindracea -81.87792 45.97071 gbif 2019-08-14 2350446287
Liatris cylindracea -90.81668 38.46749 gbif 2019-08-17 2366065927
Liatris cylindracea -87.80143 42.45080 gbif 2019-08-17 2366094074
Liatris cylindracea -87.80476 42.41891 gbif 2019-08-18 2366108212
Liatris cylindracea -87.20401 41.62109 gbif 2019-08-20 2366133525
Liatris cylindracea -81.98693 42.28380 gbif 2019-08-14 2366142035
Liatris cylindracea -81.88928 42.31200 gbif 2019-08-14 2366142290
Liatris cylindracea -89.34788 41.87327 gbif 2019-08-21 2397511006
Liatris cylindracea -89.20761 42.21227 gbif 2019-08-22 2397511044
Liatris cylindracea -89.77283 43.14011 gbif 2019-08-24 2397544031
Liatris cylindracea -81.82863 43.38344 gbif 2019-08-24 2397545195
Liatris cylindracea -80.17430 42.58551 gbif 2019-08-23 2397546505
Liatris cylindracea -88.31720 42.19544 gbif 2019-08-23 2397561898
Liatris cylindracea -88.25443 42.01464 gbif 2019-08-23 2397564596
Liatris cylindracea -81.92915 42.38429 gbif 2019-08-26 2397581937
Liatris cylindracea -81.84695 42.23364 gbif 2019-08-25 2397581984
Liatris cylindracea -81.98480 42.29205 gbif 2019-08-26 2397584078
Liatris cylindracea -82.69442 45.78553 gbif 2019-08-27 2397595498
Liatris cylindracea -89.77349 43.13907 gbif 2019-08-08 2397596192
Liatris cylindracea -91.76531 43.80910 gbif 2019-08-22 2397605137
Liatris cylindracea -95.10042 43.78117 gbif 2019-08-28 2397610615
Liatris cylindracea -87.80180 42.45135 gbif 2019-08-28 2397623166
Liatris cylindracea -89.17146 42.13536 gbif 2019-08-27 2397626885
Liatris cylindracea -89.62310 42.27654 gbif 2019-08-29 2397662767
Liatris cylindracea -81.99679 43.29185 gbif 2019-08-31 2397681829
Liatris cylindracea -82.74137 45.78097 gbif 2019-08-10 2422875270
Liatris cylindracea -89.44227 43.23465 gbif 2019-08-27 2422995129
Liatris cylindracea -86.03548 42.60680 gbif 2019-08-15 2423013181
Liatris cylindracea -81.90362 42.32986 gbif 2019-08-14 2429230412
Liatris cylindracea -81.80529 42.29334 gbif 2019-08-13 2445011465
Liatris cylindracea -81.99733 42.22599 gbif 2019-08-17 2445121474
Liatris cylindracea -91.42895 43.77074 gbif 2019-08-22 2445137307
Liatris cylindracea -80.07392 43.52519 gbif 2019-08-07 2465116216
Liatris cylindracea -88.20208 42.06030 gbif 2019-08-23 2465120000
Liatris cylindracea -81.83886 42.24473 gbif 2019-08-04 2517952699
Liatris cylindracea -80.32013 43.46092 gbif 2019-08-13 2517954637
Liatris cylindracea -81.97693 43.28063 gbif 2019-08-29 2517954794
Liatris cylindracea -81.64483 46.37250 gbif 2019-08-22 2540737415
Liatris cylindracea -89.20591 42.21140 gbif 2019-08-16 2576419875
Liatris cylindracea -89.39355 43.18280 gbif 2019-08-14 2580229150
Liatris cylindracea -89.30190 48.58295 gbif 2019-08-18 2813977735
Liatris cylindracea -81.89001 42.33353 gbif 2019-09-04 2422925996
Liatris cylindracea -88.27338 42.15473 gbif 2019-09-08 2422978756
Liatris cylindracea -89.05805 42.98933 gbif 2019-09-05 2422979604
Liatris cylindracea -82.99418 45.80779 gbif 2019-09-06 2422979686
Liatris cylindracea -83.01787 45.91057 gbif 2019-09-12 2423048695
Liatris cylindracea -81.80096 42.36401 gbif 2019-09-13 2423057406
Liatris cylindracea -82.76924 45.99909 gbif 2019-09-12 2423063072
Liatris cylindracea -91.23133 37.18044 gbif 2019-09-10 2423064456
Liatris cylindracea -90.77232 43.42788 gbif 2019-09-14 2423075440
Liatris cylindracea -90.77202 43.42764 gbif 2019-09-14 2423078113
Liatris cylindracea -82.70426 45.94042 gbif 2019-09-12 2423083061
Liatris cylindracea -83.07334 43.99958 gbif 2019-09-15 2423089176
Liatris cylindracea -90.19143 38.15793 gbif 2019-09-22 2423247267
Liatris cylindracea -81.90281 42.32196 gbif 2019-09-29 2429308076
Liatris cylindracea -91.01993 42.73446 gbif 2019-09-16 2445035832
Liatris cylindracea -81.86831 43.37468 gbif 2019-09-02 2517955229
Liatris cylindracea -81.88395 42.24337 gbif 2019-09-10 2517956003
Liatris cylindracea -88.27204 41.92578 gbif 2019-09-12 2540738310
Liatris cylindracea -90.18969 38.15720 gbif 2019-09-03 2542904732
Liatris cylindracea -90.62695 38.25671 gbif 2019-09-23 2573822203
Liatris cylindracea -88.18185 41.89795 gbif 2019-09-19 2576312278
Liatris cylindracea -81.84571 42.27222 gbif 2019-10-11 2429553598
Liatris cylindracea -91.22471 36.31743 gbif 2019-11-04 2451726900
Liatris cylindracea -87.54422 41.75960 gbif 2018-04-28 1847423249
Liatris cylindracea -89.49863 42.68726 gbif 2018-07-18 1883562246
Liatris cylindracea -81.88823 43.22347 gbif 2018-07-19 1883563104
Liatris cylindracea -91.41243 43.95158 gbif 2018-07-21 1883931517
Liatris cylindracea -81.90000 43.30244 gbif 2018-07-19 1883938526
Liatris cylindracea -91.40139 43.95866 gbif 2018-07-22 1883943054
Liatris cylindracea -83.75256 44.53482 gbif 2018-07-05 1883963011
Liatris cylindracea -91.01829 42.73412 gbif 2018-07-23 1890039481
Liatris cylindracea -80.09582 43.57053 gbif 2018-07-28 1890054025
Liatris cylindracea -88.25598 42.01284 gbif 2018-07-29 1890059243
Liatris cylindracea -87.57331 41.78203 gbif 2018-07-24 1890059904
Liatris cylindracea -87.54423 41.75961 gbif 2018-07-29 1890064325
Liatris cylindracea -80.08350 43.58716 gbif 2018-07-30 1890071904
Liatris cylindracea -81.83987 43.20576 gbif 2018-07-30 1890074211
Liatris cylindracea -79.55828 43.69791 gbif 2018-07-30 1890074974
Liatris cylindracea -83.60155 42.23280 gbif 2018-07-29 1890674100
Liatris cylindracea -90.75028 44.29741 gbif 2018-07-30 1890707606
Liatris cylindracea -90.75057 44.30033 gbif 2018-07-30 1890724619
Liatris cylindracea -88.25505 42.01478 gbif 2018-07-31 1891295584
Liatris cylindracea -83.69037 42.30337 gbif 2018-07-29 1899767042
Liatris cylindracea -91.40949 43.95753 gbif 2018-07-21 1914310744
Liatris cylindracea -84.49331 39.06933 gbif 2018-07-25 1920787697
Liatris cylindracea -89.49863 42.68726 gbif 2018-07-18 1990471361
Liatris cylindracea -80.05726 43.47835 gbif 2018-07-30 2465087376
Liatris cylindracea -79.53908 43.74333 gbif 2018-07-26 2540717475
Liatris cylindracea -87.57825 41.78325 gbif 2018-07-28 2626116351
Liatris cylindracea -80.12584 43.52970 gbif 2018-07-29 2813865345
Liatris cylindracea -89.33311 41.23892 gbif 2018-08-03 1890680138
Liatris cylindracea -87.80849 42.42209 gbif 2018-08-04 1890688568
Liatris cylindracea -81.90017 42.38337 gbif 2018-08-05 1890698827
Liatris cylindracea -89.05772 42.98993 gbif 2018-08-05 1890710637
Liatris cylindracea -77.84414 44.03936 gbif 2018-08-07 1890712871
Liatris cylindracea -89.77215 43.13979 gbif 2018-08-03 1890715967
Liatris cylindracea -87.87305 41.77304 gbif 2018-08-11 1891310532
Liatris cylindracea -89.38538 41.92383 gbif 2018-08-08 1891325032
Liatris cylindracea -81.84592 42.30069 gbif 2018-08-12 1891325190
Liatris cylindracea -87.57331 41.78203 gbif 2018-08-11 1891327188
Liatris cylindracea -88.31518 42.19575 gbif 2018-08-09 1891334163
Liatris cylindracea -88.31226 42.14957 gbif 2018-08-17 1898797804
Liatris cylindracea -87.80472 42.41842 gbif 2018-08-17 1898800066
Liatris cylindracea -87.80514 42.41889 gbif 2018-08-17 1898800190
Liatris cylindracea -87.80269 42.45021 gbif 2018-08-18 1898810701
Liatris cylindracea -89.05806 42.98947 gbif 2018-08-18 1898834065
Liatris cylindracea -89.77281 43.14006 gbif 2018-08-14 1898834091
Liatris cylindracea -92.65693 44.52560 gbif 2018-08-21 1898843056
Liatris cylindracea -81.98080 42.21883 gbif 2018-08-14 1901034927
Liatris cylindracea -89.05834 42.99005 gbif 2018-08-25 1901036256
Liatris cylindracea -87.80143 42.45066 gbif 2018-08-03 1914194648
Liatris cylindracea -87.78765 41.86717 gbif 2018-08-14 1914195357
Liatris cylindracea -89.36965 41.90080 gbif 2018-08-18 1914195583
Liatris cylindracea -87.80263 42.45006 gbif 2018-08-19 1914195758
Liatris cylindracea -87.57331 41.78203 gbif 2018-08-18 1914196356
Liatris cylindracea -86.04611 42.76894 gbif 2018-08-18 1914348645
Liatris cylindracea -91.07192 37.09674 gbif 2018-08-31 1914352037
Liatris cylindracea -87.80474 42.41722 gbif 2018-08-22 1948753365
Liatris cylindracea -79.50747 43.75722 gbif 2018-08-19 2005266956
Liatris cylindracea -91.81193 38.14330 gbif 2018-08-13 2006091296
Liatris cylindracea -81.88098 43.20848 gbif 2018-08-19 2294401867
Liatris cylindracea -83.07481 43.99593 gbif 2018-08-05 2397448175
Liatris cylindracea -82.69326 45.73200 gbif 2018-08-07 2422915938
Liatris cylindracea -82.66572 45.77302 gbif 2018-08-07 2422916195
Liatris cylindracea -80.15631 43.44793 gbif 2018-08-07 2465087402
Liatris cylindracea -83.45101 38.88592 gbif 2018-08-04 2529316034
Liatris cylindracea -81.88460 43.29642 gbif 2018-08-04 2540717651
Liatris cylindracea -80.46163 42.73060 gbif 2018-08-01 2557766806
Liatris cylindracea -89.06015 42.99205 gbif 2018-09-08 1913164558
Liatris cylindracea -89.05841 42.98910 gbif 2018-09-08 1913165897
Liatris cylindracea -81.98628 42.23625 gbif 2018-09-08 1913178556
Liatris cylindracea -88.25508 42.01475 gbif 2018-09-16 1913685389
Liatris cylindracea -90.54347 38.20224 gbif 2018-10-01 1914345018
Liatris cylindracea -87.80004 42.45454 gbif 2018-10-06 1920805907
Liatris cylindracea -80.59191 42.70125 gbif 2018-10-13 1931495094
Liatris cylindracea -87.74641 40.21564 gbif 2018-10-14 1931503386
Liatris cylindracea -87.74547 40.21462 gbif 2018-10-14 1931503491
Liatris cylindracea -87.93044 42.57168 gbif 2018-10-13 1932309128
Liatris cylindracea -89.74419 44.27951 gbif 2017-07-09 1572386698
Liatris cylindracea -79.57448 43.67619 gbif 2017-07-12 1586098036
Liatris cylindracea -88.48202 42.87147 gbif 2017-07-27 1586107794
Liatris cylindracea -89.05739 42.99125 gbif 2017-07-14 1640041333
Liatris cylindracea -89.05973 42.99013 gbif 2017-07-23 1640081143
Liatris cylindracea -89.80954 43.10593 gbif 2017-07-30 1640088544
Liatris cylindracea -89.39449 43.18241 gbif 2017-07-30 1668883155
Liatris cylindracea -88.36231 42.09077 gbif 2017-07-03 1836687918
Liatris cylindracea -89.01868 42.32850 gbif 2017-08-06 1586144059
Liatris cylindracea -90.51111 41.43678 gbif 2017-08-07 1586148475
Liatris cylindracea -84.65049 44.65616 gbif 2017-08-04 1586925946
Liatris cylindracea -87.54650 41.75860 gbif 2017-08-12 1586934321
Liatris cylindracea -90.91875 43.03481 gbif 2017-08-12 1586936620
Liatris cylindracea -90.59550 43.17681 gbif 2017-08-12 1586938713
Liatris cylindracea -90.59428 43.17674 gbif 2017-08-12 1586939642
Liatris cylindracea -89.05808 42.98948 gbif 2017-08-14 1586947741
Liatris cylindracea -87.80483 42.41762 gbif 2017-08-19 1621782362
Liatris cylindracea -89.06039 42.99199 gbif 2017-08-18 1621782908
Liatris cylindracea -90.80156 44.32064 gbif 2017-08-20 1621801032
Liatris cylindracea -77.98699 44.08402 gbif 2017-08-23 1621802198
Liatris cylindracea -77.97258 44.18646 gbif 2017-08-23 1621802202
Liatris cylindracea -77.84793 44.11483 gbif 2017-08-23 1621802338
Liatris cylindracea -89.12631 40.65682 gbif 2017-08-20 1640049023
Liatris cylindracea -89.77192 43.13977 gbif 2017-08-25 1640090810
Liatris cylindracea -89.51541 43.02454 gbif 2017-08-06 1640092507
Liatris cylindracea -89.12686 40.65733 gbif 2017-08-23 1640099964
Liatris cylindracea -90.59409 43.17872 gbif 2017-08-12 1640104348
Liatris cylindracea -90.92108 43.03458 gbif 2017-08-12 1640119043
Liatris cylindracea -89.05958 42.99066 gbif 2017-08-26 1640121164
Liatris cylindracea -89.05909 42.99085 gbif 2017-08-26 1668826441
Liatris cylindracea -81.98498 43.23406 gbif 2017-08-28 1668843264
Liatris cylindracea -89.05973 42.99013 gbif 2017-08-05 1677281046
Liatris cylindracea -87.80461 42.41739 gbif 2017-08-13 1677326118
Liatris cylindracea -89.18842 42.93090 gbif 2017-08-05 1703207537
Liatris cylindracea -86.30509 35.85191 gbif 2017-08-25 1914197983
Liatris cylindracea -77.84368 44.14095 gbif 2017-08-10 1993753945
Liatris cylindracea -83.41716 38.97641 gbif 2017-08-19 2236221183
Liatris cylindracea -83.46504 38.80499 gbif 2017-08-19 2236221520
Liatris cylindracea -89.05756 42.99076 gbif 2017-08-05 2445094349
Liatris cylindracea -92.16922 35.96823 gbif 2017-09-03 1640089470
Liatris cylindracea -77.86183 44.01423 gbif 2017-09-05 1640099376
Liatris cylindracea -89.49340 42.69841 gbif 2017-09-06 1640104962
Liatris cylindracea -89.71798 43.13156 gbif 2017-09-08 1640120915
Liatris cylindracea -77.88411 44.05421 gbif 2017-09-05 1699378358
Liatris cylindracea -89.39243 43.18356 gbif 2017-09-19 1831149516
Liatris cylindracea -87.57608 41.71460 gbif 2016-07-26 1291154225
Liatris cylindracea -87.64774 41.84284 gbif 2016-07-28 1291154355
Liatris cylindracea -89.05757 42.99036 gbif 2016-07-30 1291159304
Liatris cylindracea -85.71997 43.44973 gbif 2016-07-23 1322997904
Liatris cylindracea -87.99824 42.52563 gbif 2016-07-30 1453168486
Liatris cylindracea -87.90689 42.48217 gbif 2016-07-30 1453168650
Liatris cylindracea -83.36147 45.96419 gbif 2016-07-24 2012952611
Liatris cylindracea -82.73759 45.61861 gbif 2016-07-24 2012952614
Liatris cylindracea -77.86305 44.07157 gbif 2016-07-31 2563557365
Liatris cylindracea -87.82715 42.36894 gbif 2016-08-06 1291670420
Liatris cylindracea -89.33770 41.90610 gbif 2016-08-08 1291677710
Liatris cylindracea -89.06095 42.99140 gbif 2016-08-20 1315052190
Liatris cylindracea -80.06046 42.54193 gbif 2016-08-11 2563549536
Liatris cylindracea -89.84055 40.14877 gbif 2016-09-05 1306585813
Liatris cylindracea -95.25413 45.44381 gbif 2015-07-16 1805364817
Liatris cylindracea -92.06722 44.47028 gbif 2015-08-10 1135214267
Liatris cylindracea -83.45041 38.90135 gbif 2015-08-19 1143523705
Liatris cylindracea -77.89793 44.06960 gbif 2015-08-09 1453190881
Liatris cylindracea -81.86100 42.20755 gbif 2015-08-22 1978431574
Liatris cylindracea -86.20510 44.16520 gbif 2015-08-13 1989262155
Liatris cylindracea -82.03146 45.55013 gbif 2015-08-29 1990527406
Liatris cylindracea -90.62555 38.25742 gbif 2015-09-23 1262228557
Liatris cylindracea -89.06000 42.99222 gbif 2015-09-02 1455566941
Liatris cylindracea -86.30920 35.84616 gbif 2015-09-08 1455566952
Liatris cylindracea -92.12878 35.95507 gbif 2015-09-07 1455566954
Liatris cylindracea -88.05258 35.46717 gbif 2015-09-11 2609176513
Liatris cylindracea -89.39301 43.18290 gbif 2014-03-30 1455597299
Liatris cylindracea -89.39100 43.18405 gbif 2014-03-30 1455597300
Liatris cylindracea -86.04586 42.60770 gbif 2014-07-29 1988722401
Liatris cylindracea -90.67639 44.19250 gbif 2014-08-02 1024203052
Liatris cylindracea -88.20065 42.06067 gbif 2014-08-15 1024209902
Liatris cylindracea -84.64898 44.62078 gbif 2014-08-26 1024218875
Liatris cylindracea -90.82306 38.46583 gbif 2014-08-19 1258065292
Liatris cylindracea -80.11498 42.48124 gbif 2014-08-26 1703218792
Liatris cylindracea -87.02639 33.05893 gbif 2014-08-17 2005285277
Liatris cylindracea -82.95888 45.98204 gbif 2014-09-02 2012954819
Liatris cylindracea -77.81174 44.19346 gbif 2014-09-05 2563557752
Liatris cylindracea -90.52833 37.97622 gbif 2013-07-26 1257980954
Liatris cylindracea -83.09683 45.93570 gbif 2013-07-29 1931528898
Liatris cylindracea -92.10115 44.39779 gbif 2013-08-23 1990528787
Liatris cylindracea -91.70500 36.65830 gbif 2013-09-08 1257994504
Liatris cylindracea -80.03283 42.56807 gbif 2013-09-14 2294728102
Liatris cylindracea -88.00617 35.54113 gbif 2013-09-10 2609336706
Liatris cylindracea -81.98543 42.22754 gbif 2013-09-08 899942247
Liatris cylindracea -91.36310 37.18640 gbif 2013-10-23 1258050633
Liatris cylindracea -92.44970 36.69610 gbif 2013-10-22 1258050652
Liatris cylindracea -82.71642 45.67176 gbif 2012-07-14 1890712775
Liatris cylindracea -81.93289 43.20522 gbif 2012-08-19 1802599909
Liatris cylindracea -88.00613 35.54085 gbif 2012-08-18 2540715851
Liatris cylindracea -90.74426 44.29697 gbif 2012-08-04 891781755
Liatris cylindracea -90.54080 38.20420 gbif 2012-09-24 1257883005
Liatris cylindracea -89.77221 43.13984 gbif 2011-08-14 1455597190
Liatris cylindracea -80.07017 42.55773 gbif 2011-08-20 1990508689
Liatris cylindracea -87.39828 41.62251 gbif 2010-07-17 891748737
Liatris cylindracea -80.19311 42.43746 gbif 2010-08-02 1836586916
Liatris cylindracea -91.03444 38.31917 gbif 2009-08-21 2268820435
Liatris cylindracea -82.64870 45.74304 gbif 2009-08-26 2366085536
Liatris cylindracea -79.51408 43.74617 gbif 2009-10-20 2265807122
Liatris cylindracea -83.11034 45.93447 gbif 2007-07-19 2005339463
Liatris cylindracea -87.74690 40.21827 gbif 2007-07-09 2234780444
Liatris cylindracea -80.05904 42.46448 gbif 2007-09-08 1677316926
Liatris cylindracea -92.72222 36.80083 gbif 2006-09-08 1260453490
Liatris cylindracea -90.70861 38.67861 gbif 2004-08-30 1261008488
Liatris cylindracea -87.80857 42.42208 gbif 2003-08-17 1453373091
Liatris cylindracea -81.96000 43.21000 gbif 2002-08-14 1456194679
Liatris cylindracea -81.96000 43.21000 gbif 2002-08-14 1899312387
Liatris cylindracea -81.96000 43.21000 gbif 2002-08-14 2516828214
Liatris cylindracea -93.66378 45.45102 gbif 2001-08-28 2265430701
Liatris cylindracea -90.54083 38.20388 gbif 2001-09-20 1258814089
Liatris cylindracea -93.72596 45.46555 gbif 1998-08-14 2265349667
Liatris cylindracea -91.42269 43.75640 gbif 1998-08-08 2265430141
Liatris cylindracea -91.42269 43.75640 gbif 1997-08-07 2265443698
Liatris cylindracea -90.50000 37.78333 gbif 1996-09-23 1258330970
Liatris cylindracea -91.03333 37.16667 gbif 1995-08-01 1261880523
Liatris cylindracea -92.01953 43.76751 gbif 1994-07-29 2265349375
Liatris cylindracea -92.08925 44.09956 gbif 1994-07-26 2265350368
Liatris cylindracea -89.87222 40.69361 gbif 1994-08-18 415908558
Liatris cylindracea -91.70194 36.77805 gbif 1993-08-15 1260978526
Liatris cylindracea -90.58333 37.96667 gbif 1993-08-11 1261296712
Liatris cylindracea -104.96147 39.73221 gbif 1993-08-17 2242400797
Liatris cylindracea -92.74229 44.63663 gbif 1993-08-04 2265349284
Liatris cylindracea -92.74229 44.63663 gbif 1993-08-02 2265349389
Liatris cylindracea -87.10000 33.00000 gbif 1992-07-19 56603321
Liatris cylindracea -90.67444 38.12500 gbif 1992-08-15 2268932005
Liatris cylindracea -90.67444 38.12500 gbif 1992-08-15 2464647947
Liatris cylindracea -91.11611 38.18805 gbif 1992-09-24 1261002493
Liatris cylindracea -84.55462 44.74868 gbif 1992-09-02 1988989334
Liatris cylindracea -84.57184 44.69566 gbif 1992-09-01 1989228152
Liatris cylindracea -90.67500 38.12500 gbif 1992-10-05 2268929729
Liatris cylindracea -86.03300 42.60900 gbif 1991-08-08 233099622
Liatris cylindracea -90.70000 38.16666 gbif 1990-08-28 1261312058
Liatris cylindracea -90.62055 38.25888 gbif 1990-09-02 1261003426
Liatris cylindracea -93.19351 45.39232 gbif 1989-08-15 2265349359
Liatris cylindracea -92.54074 44.49194 gbif 1988-06-15 2265349858
Liatris cylindracea -93.53888 36.51083 gbif 1988-07-25 1260965018
Liatris cylindracea -92.81944 44.80028 gbif 1988-07-26 2265349442
Liatris cylindracea -84.54181 44.45902 gbif 1988-08-22 1989020896
Liatris cylindracea -92.91666 36.68333 gbif 1987-09-19 1261002967
Liatris cylindracea -92.79092 44.84076 gbif 1987-09-30 2265349996
Liatris cylindracea -92.10861 38.41472 gbif 1986-07-18 1260964907
Liatris cylindracea -91.99041 44.17096 gbif 1986-07-31 2265349957
Liatris cylindracea -92.54583 44.52361 gbif 1986-08-11 2265349388
Liatris cylindracea -92.77505 44.85476 gbif 1985-06-02 2265349779
Liatris cylindracea -92.54166 38.11666 gbif 1985-07-25 1261007084
Liatris cylindracea -80.83330 40.00000 gbif 1985-08-12 1228196013
Liatris cylindracea -83.43154 44.44687 gbif 1985-08-01 1988824682
Liatris cylindracea -87.00445 33.00889 gbif 1985-08-01 56603146
Liatris cylindracea -95.24502 45.44867 gbif 1982-07-29 2265349246
Liatris cylindracea -92.52110 44.56198 gbif 1982-07-31 2265349484
Liatris cylindracea -83.51282 44.54746 gbif 1982-08-01 1989037167
Liatris cylindracea -83.41114 44.44674 gbif 1982-08-07 1989165286
Liatris cylindracea -92.04933 44.04203 gbif 1982-08-16 2265349631
Liatris cylindracea -82.75000 46.78333 gbif 1981-08-19 1804548575
Liatris cylindracea -81.86667 46.00000 gbif 1981-08-15 1804560873
Liatris cylindracea -91.99041 44.17096 gbif 1980-08-15 2265350478
Liatris cylindracea -96.44879 47.73852 gbif 1979-08-15 2265350144
Liatris cylindracea -91.42281 43.76774 gbif 1978-08-14 2265349218
Liatris cylindracea -91.42233 43.62329 gbif 1978-08-15 2265349369
Liatris cylindracea -81.83333 43.26667 gbif 1978-09-01 1804428271
Liatris cylindracea -81.83333 43.26667 gbif 1978-09-01 1804574090
Liatris cylindracea -81.83333 43.26667 gbif 1978-09-01 1804600129
Liatris cylindracea -92.88527 36.69166 gbif 1977-07-15 1260946785
Liatris cylindracea -84.65981 45.36937 gbif 1977-07-12 1988815636
Liatris cylindracea -82.82940 45.80870 gbif 1976-07-16 1839334600
Liatris cylindracea -80.83492 45.27013 gbif 1976-08-05 1839334611
Liatris cylindracea -81.86621 42.28171 gbif 1976-09-05 1839334612
Liatris cylindracea -93.13361 35.27833 gbif 1972-10-14 924820107
Liatris cylindracea -93.13361 35.27833 gbif 1972-10-14 924820112
Liatris cylindracea -85.94579 42.12044 gbif 1970-08-29 1989329755
Liatris cylindracea -90.67583 38.50361 gbif 1969-08-01 1260932129
Liatris cylindracea -79.42300 43.85700 gbif 1968-08-25 2250594823
Liatris cylindracea -79.42000 43.86000 gbif 1968-08-25 2304990669
Liatris cylindracea -79.42000 43.86000 gbif 1968-08-25 2309958506
Liatris cylindracea -84.36783 42.05598 gbif 1967-08-15 1988659685
Liatris cylindracea -91.36237 43.62304 gbif 1966-07-29 2265350030
Liatris cylindracea -87.59237 41.56148 gbif 1963-08-15 2242401085
Liatris cylindracea -91.55129 43.63067 gbif 1963-08-21 2265349140
Liatris cylindracea -80.60000 42.90000 gbif 1963-09-19 2250594358
Liatris cylindracea -80.60000 42.90000 gbif 1963-09-19 2305014115
Liatris cylindracea -80.60000 42.90000 gbif 1963-09-19 2309823689
Liatris cylindracea -81.54317 44.34863 gbif 1962-07-30 1638340768
Liatris cylindracea -81.54317 44.34863 gbif 1962-07-30 1839334604
Liatris cylindracea -92.49944 36.92528 gbif 1961-08-02 1260987544
Liatris cylindracea -92.52110 44.56198 gbif 1960-08-21 2265350298
Liatris cylindracea -82.38690 42.99300 gbif 1957-06-09 2250594555
Liatris cylindracea -82.39000 42.99000 gbif 1957-06-09 2305853966
Liatris cylindracea -82.39000 42.99000 gbif 1957-06-09 2308553621
Liatris cylindracea -83.69666 44.44594 gbif 1957-08-01 1989178809
Liatris cylindracea -85.80775 44.16614 gbif 1957-08-18 1989380604
Liatris cylindracea -91.22528 39.55194 gbif 1956-10-06 1260987611
Liatris cylindracea -85.81486 42.13144 gbif 1955-09-07 1988830248
Liatris cylindracea -84.31931 44.41350 gbif 1954-07-12 1988708823
Liatris cylindracea -84.40077 44.50095 gbif 1953-07-25 1989311280
Liatris cylindracea -84.70668 44.63540 gbif 1953-07-24 1989351547
Liatris cylindracea -84.12964 44.57944 gbif 1953-09-27 1988687346
Liatris cylindracea -83.69494 44.34345 gbif 1952-07-20 1989524892
Liatris cylindracea -84.03907 44.57348 gbif 1951-07-28 1989202098
Liatris cylindracea -84.28180 44.61553 gbif 1951-07-31 1989288026
Liatris cylindracea -83.58146 44.05117 gbif 1951-08-11 1988784418
Liatris cylindracea -83.41105 44.42582 gbif 1951-08-11 1989049524
Liatris cylindracea -83.20592 43.97216 gbif 1951-08-09 1989177001
Liatris cylindracea -84.19955 42.29950 gbif 1951-09-15 1989290190
Liatris cylindracea -86.12800 43.35985 gbif 1949-07-30 1989187265
Liatris cylindracea -85.62947 43.02106 gbif 1948-08-17 1988599605
Liatris cylindracea -81.75590 43.31246 gbif 1945-09-04 465928982
Liatris cylindracea -91.39316 43.95302 gbif 1943-08-25 2265350350
Liatris cylindracea -79.42000 43.70000 gbif 1942-08-14 1988065648
Liatris cylindracea -85.88245 43.90066 gbif 1941-08-03 1989207965
Liatris cylindracea -84.00525 42.45669 gbif 1941-08-16 1989394120
Liatris cylindracea -85.97790 42.54850 gbif 1939-07-22 1989545608
Liatris cylindracea -85.97790 42.54850 gbif 1939-08-07 1989365814
Liatris cylindracea -80.24910 42.55910 gbif 1938-08-25 1839334613
Liatris cylindracea -92.65777 37.08861 gbif 1937-07-28 1260970503
Liatris cylindracea -92.21972 38.43194 gbif 1937-08-18 1260971298
Liatris cylindracea -91.99472 37.91638 gbif 1937-08-26 1260971511
Liatris cylindracea -92.17861 37.89166 gbif 1937-08-26 1260971530
Liatris cylindracea -91.79694 37.74083 gbif 1937-08-27 1260971557
Liatris cylindracea -91.00083 39.21722 gbif 1937-09-08 1260971843
Liatris cylindracea -91.91861 38.80444 gbif 1937-09-10 1260971923
Liatris cylindracea -91.31256 43.54390 gbif 1936-08-20 2265349387
Liatris cylindracea -79.46667 43.63333 gbif 1934-08-20 1899286238
Liatris cylindracea -79.52797 43.66711 gbif 1934-08-20 465928980
Liatris cylindracea -82.08000 42.75000 gbif 1934-09-02 1988141918
Liatris cylindracea -81.50000 43.67000 gbif 1934-09-03 1988211485
Liatris cylindracea -90.44777 38.56833 gbif 1931-08-28 1260931882
Liatris cylindracea -79.48333 43.63333 gbif 1929-09-15 1899285105
Liatris cylindracea -91.97750 37.92611 gbif 1928-09-02 2268815320
Liatris cylindracea -91.97750 37.92611 gbif 1927-09-09 2268815131
Liatris cylindracea -83.13098 42.67714 gbif 1925-08-26 1988722234
Liatris cylindracea -83.72684 42.29218 gbif 1923-09-01 1989085003
Liatris cylindracea -91.57530 36.98690 gbif 1922-08-04 2268816553
Liatris cylindracea -79.46667 43.63333 gbif 1920-08-07 1899369426
Liatris cylindracea -84.61130 45.37999 gbif 1920-08-21 1989176004
Liatris cylindracea -91.57528 36.98694 gbif 1920-10-05 1260957420
Liatris cylindracea -91.31853 43.62916 gbif 1919-08-01 2265351238
Liatris cylindracea -83.61299 42.24115 gbif 1917-09-03 1989001083
Liatris cylindracea -91.49750 40.04277 gbif 1915-09-06 1260939070
Liatris cylindracea -91.97750 37.92611 gbif 1913-10-16 2268815553
Liatris cylindracea -91.57530 36.98690 gbif 1910-08-05 2268813661
Liatris cylindracea -93.16646 44.90516 gbif 1909-08-07 2265349893
Liatris cylindracea -91.57530 36.98690 gbif 1908-09-10 2268814511
Liatris cylindracea -91.57530 36.98690 gbif 1905-10-08 2268813695
Liatris cylindracea -91.28287 43.50806 gbif 1899-07-29 2265350102
Liatris cylindracea -96.33441 30.62798 gbif 1897-01-01 1930308135
Liatris cylindracea -93.75416 36.54805 gbif 1896-09-28 1260926595
Liatris cylindracea -90.67583 38.50361 gbif 1896-09-10 2268814493
Liatris cylindracea -91.15305 41.46611 gbif 1894-08-18 415908555
Liatris cylindracea -83.72633 42.27087 gbif 1892-08-12 1989425627
Liatris cylindracea -92.26684 44.44941 gbif 1892-08-20 2265350910
Liatris cylindracea -91.63931 44.04996 gbif 1888-09-13 2265350426
Liatris cylindracea -92.36995 44.41148 gbif 1883-08-14 2265350636
Liatris cylindracea -90.67583 38.50361 gbif 1883-08-25 2268815378
Liatris cylindracea -84.83763 43.09254 gbif 1880-01-01 1988606890
Liatris cylindracea -88.99055 40.51416 gbif 1879-08-30 415908557
Liatris cylindracea -92.98853 42.11249 gbif 1878-08-01 1928015915
Liatris cylindracea -84.83763 43.09254 gbif 1878-08-01 1989045372
Liatris cylindracea -93.26384 44.97996 gbif 1878-08-14 2265350145
Liatris cylindracea -90.29208 38.64010 gbif 1875-08-30 1260941240
Liatris cylindracea -83.30227 42.51821 gbif 1874-01-01 1988958465
Liatris cylindracea -83.72633 42.27087 gbif 1861-07-31 1988611173
Liatris cylindracea -83.72633 42.27087 gbif 1861-08-23 1989097124
Liatris cylindracea -79.00796 43.09099 gbif 1852-08-16 1928848145
Liatris cylindracea -91.05139 39.44889 gbif 1850-01-01 1258836178
Liatris cylindracea -80.76667 34.76667 gbif NA 1144624572
Liatris cylindracea -79.80561 43.89389 gbif NA 1899385146
Liatris cylindracea -83.72633 42.27087 gbif NA 1989409051
Liatris cylindracea -92.83619 45.42608 gbif NA 2265349156
Liatris cylindracea -95.38976 45.65024 gbif NA 2265349165
Liatris cylindracea -91.74803 43.99779 gbif NA 2265349170
Liatris cylindracea -92.72282 44.40976 gbif NA 2265349268
Liatris cylindracea -95.38976 45.65024 gbif NA 2265349376
Liatris cylindracea -92.66120 44.28943 gbif NA 2265349481
Liatris cylindracea -93.16160 44.45830 gbif NA 2265349647
Liatris cylindracea -93.28313 45.54365 gbif NA 2265349799
Liatris cylindracea -92.61081 44.49954 gbif NA 2265349807
Liatris cylindracea -93.49197 44.97475 gbif NA 2265349841
Liatris cylindracea -93.28313 45.54365 gbif NA 2265349855
Liatris cylindracea -92.99108 45.50656 gbif NA 2265349865
Liatris cylindracea -95.38976 45.65024 gbif NA 2265349960
Liatris cylindracea -94.94038 45.19606 gbif NA 2265350037
Liatris cylindracea -94.24059 46.34845 gbif NA 2265350074
Liatris cylindracea -93.26384 44.97996 gbif NA 2265350154
Liatris cylindracea -93.59000 44.02000 gbif NA 2265350197
Liatris cylindracea -92.72259 44.40985 gbif NA 2265350334
Liatris cylindracea -93.16693 45.08003 gbif NA 2265350357
Liatris cylindracea -93.26384 44.97996 gbif NA 2265350396
Liatris cylindracea -93.47000 45.00000 gbif NA 2265350467
Liatris cylindracea -95.37639 45.65025 gbif NA 2265442926
Liatris cylindracea -94.20071 46.36871 gbif NA 2265443303
Liatris cylindracea -83.53028 38.81083 gbif NA 91660785
Liatris cylindracea -87.03444 33.05861 gbif NA 91660797
Liatris cylindracea -92.83619 45.42608 idigbio NA 021b6c1d-9b22-45e4-9cea-a26821e65910
Liatris cylindracea -91.31853 43.62916 idigbio 1919-08-01 039f4bde-2ab1-4561-a9af-1c850a42aab9
Liatris cylindracea -86.00100 38.86488 idigbio 1958-09-07 0923ee62-4829-4879-a6ae-d8cfec88e331
Liatris cylindracea -91.42281 43.76774 idigbio 1978-08-14 09848e55-039e-49b8-91f5-7c3cded80e1c
Liatris cylindracea -90.57678 41.52914 idigbio NA 0bab91db-f7db-4115-aebf-18bd0533cf25
Liatris cylindracea -91.42233 43.62329 idigbio 1978-08-15 0bbf443b-785b-4dde-8421-77c60fdfacac
Liatris cylindracea -93.26384 44.97996 idigbio 1878-08-14 0fafba30-3764-4fd9-aed5-f3e98531721c
Liatris cylindracea -91.99041 44.17096 idigbio 1980-08-15 103acddc-33b1-4ea4-a816-47eb335aea70
Liatris cylindracea -92.54583 44.52361 idigbio 1986-08-11 134b3de9-aa79-4e23-8e1e-1a8dd0149f68
Liatris cylindracea -83.72633 42.27087 idigbio 1861-08-23 13b34217-ea49-468e-8bfe-ccb255446415
Liatris cylindracea -92.08925 44.09956 idigbio 1994-07-26 153516e5-dcba-4e68-b6cb-2e087244260b
Liatris cylindracea -89.46750 31.12278 idigbio 2011-10-31 15c9ead4-210a-41c1-921b-ba8aeb23fd81
Liatris cylindracea -84.65981 45.36937 idigbio 1977-07-12 1e31791c-986b-464b-b3f0-ddab33d5abae
Liatris cylindracea -92.04933 44.04203 idigbio 1982-08-16 20960f79-70e2-4c12-b416-f35af7198cfd
Liatris cylindracea -93.26384 44.97996 idigbio NA 239b13e3-4551-4a79-8688-462cdf4c5300
Liatris cylindracea -93.19351 45.39232 idigbio 1989-08-15 2a1eb429-d266-4d2c-a175-97c04cd19ea8
Liatris cylindracea -81.83333 43.26667 idigbio 1978-09-01 2a826fb4-47bf-401f-bcf6-b0c42f3f4654
Liatris cylindracea -83.72633 42.27087 idigbio 1892-08-12 2a9abca5-0b70-479f-a360-af1dbb027b01
Liatris cylindracea -93.16160 44.45830 idigbio NA 2a9d5de3-06a6-4f91-be67-40344700ca34
Liatris cylindracea -91.55129 43.63067 idigbio 1963-08-21 2c2f293c-b418-4a8f-9edb-f8fe14ea8a3e
Liatris cylindracea -81.96000 43.21000 idigbio 2002-08-14 2cdf4ab4-a796-4671-bd37-0d7e6484a1d0
Liatris cylindracea -91.08437 41.33353 idigbio 1911-07-14 323e9ef9-9f91-47ef-b78f-c619162c7367
Liatris cylindracea -92.52110 44.56198 idigbio 1960-08-21 32cd4c8c-8b98-4b89-879f-bb3e620cc57c
Liatris cylindracea -90.80709 41.45920 idigbio NA 35857481-d738-443d-816e-c3c2f22cb7a2
Liatris cylindracea -92.26684 44.44941 idigbio 1892-08-20 3d9f82c3-26b3-4e8d-96a5-36ca3c3e693d
Liatris cylindracea -81.83333 43.26667 idigbio 1978-09-01 3fd73215-55e8-45e9-863f-b3e5fc6a2c23
Liatris cylindracea -95.38976 45.65024 idigbio NA 42e283a0-0ffe-441f-ac4d-e86352a4fcca
Liatris cylindracea -88.66035 33.41494 idigbio 1992-05-27 4658101d-9236-40f6-98b2-2d8dc877abd2
Liatris cylindracea -93.26384 44.97996 idigbio NA 49d68db4-b8ce-41a6-bba8-e18f0feb2225
Liatris cylindracea -92.99108 45.50656 idigbio NA 49dcdb38-7af9-48ca-9d95-614b0e75471a
Liatris cylindracea -93.16646 44.90516 idigbio 1909-08-07 4b18242d-e06a-48b2-8646-208e536b65a5
Liatris cylindracea -86.03300 42.60900 idigbio 1991-08-08 5260420f-c130-4938-90ac-c8cffb753ba4
Liatris cylindracea -91.74803 43.99779 idigbio NA 5349f5b8-36c9-4cb3-aa1b-b9c2f0133453
Liatris cylindracea -94.94038 45.19606 idigbio NA 5364a7e5-d11d-4e9e-94a1-3a95b0eb2f80
Liatris cylindracea -90.57678 41.52914 idigbio NA 5805d2a2-ab42-4406-987e-30a20880f3c2
Liatris cylindracea -95.38976 45.65024 idigbio NA 5e2d6af4-9b80-41e3-9461-1a1f5e011bd2
Liatris cylindracea -89.64108 41.35768 idigbio 2012-08-24 61449c5e-3214-4bbe-9545-2aa02ad6c835
Liatris cylindracea -81.86667 46.00000 idigbio 1981-08-15 65596114-c374-44cb-a0f5-c785a1fd5fae
Liatris cylindracea -92.54074 44.49194 idigbio 1988-06-15 673d05ee-a10f-40d5-907e-d75b9f1de5eb
Liatris cylindracea -104.96147 39.73221 idigbio 1993-08-17 6740f27b-3023-4b34-897b-c64b82a65f7b
Liatris cylindracea -81.83333 43.26667 idigbio 1978-09-01 6bdf0388-5175-4e0a-9198-e40fd2904e92
Liatris cylindracea -83.30227 42.51821 idigbio NA 6c162c78-3483-40fb-86a7-83bac5adac7a
Liatris cylindracea -95.38976 45.65024 idigbio NA 6c70a145-7efb-419f-a0ea-02503124143a
Liatris cylindracea -92.61081 44.49954 idigbio NA 6fb0eac5-7072-40ad-a4d1-adf9a0d98338
Liatris cylindracea -92.79092 44.84076 idigbio 1987-09-30 741ead93-4cdc-42d6-aeea-f3a3ea992436
Liatris cylindracea -93.47000 45.00000 idigbio NA 765464bd-e157-4595-9cd7-6f46cbf76c86
Liatris cylindracea -83.61299 42.24115 idigbio 1917-09-03 786e298b-4fff-4379-949d-7217a3cabf48
Liatris cylindracea -91.31256 43.54390 idigbio 1936-08-20 7984f7d7-bb49-462d-9b73-d00f0324e56c
Liatris cylindracea -83.20592 43.97216 idigbio 1951-08-09 7d71d2a3-a765-457b-b52d-987b101577bc
Liatris cylindracea -88.30687 42.23804 idigbio 2011-08-03 815f35ad-4b0f-49c4-a99e-dcfa11be1585
Liatris cylindracea -92.74229 44.63663 idigbio 1993-08-02 8853b35e-736e-4b77-8cb6-569cfad0a4ca
Liatris cylindracea -91.42269 43.75640 idigbio 1997-08-07 8c48b775-d4a2-4332-b1d9-68bfdcc72dc7
Liatris cylindracea -90.55552 41.54769 idigbio 1889-07-07 8f42912a-cf30-4afc-830f-5598b79f22ee
Liatris cylindracea -86.01568 42.52116 idigbio 1939-07-22 8fac66f4-91f3-421a-82da-6d425901b846
Liatris cylindracea -91.08168 41.32893 idigbio 1910-10-01 9070fedd-e94e-4339-bd7b-cde0a538dc6d
Liatris cylindracea -84.40077 44.50095 idigbio 1953-07-25 92441d82-8af3-4635-86a5-8fe986dc5eed
Liatris cylindracea -71.09533 42.52565 idigbio 1880-08-21 9375e4f5-212e-4f1c-bfcc-cc9a0b45096a
Liatris cylindracea -89.61468 41.35802 idigbio 2010-07-07 96c2fda8-8a2e-403f-82f9-c5dadeb20c1f
Liatris cylindracea -93.13361 35.27833 idigbio 1972-10-14 9812a3fb-fe31-4568-99a4-726dad5a9164
Liatris cylindracea -86.01568 42.52116 idigbio 1939-08-07 9b953e07-596d-4fdd-b405-ffe0403e6431
Liatris cylindracea -89.46750 31.12278 idigbio 2011-10-31 9bce918d-84b7-4453-b6eb-3c176453a11b
Liatris cylindracea -93.66378 45.45102 idigbio 2001-08-28 9cd07271-ccfa-46fa-9813-daa2ff223aa7
Liatris cylindracea -94.24059 46.34845 idigbio NA a39cb742-426d-4f51-af35-e3199af1b439
Liatris cylindracea -88.09183 41.81472 idigbio 2012-09-14 a6bdfa20-fdea-4e7d-8e71-6c91db2cc44f
Liatris cylindracea -90.95013 40.23004 idigbio NA a7b08fba-3408-4fb1-ad95-370130bb0e52
Liatris cylindracea -93.49197 44.97475 idigbio NA a80e154c-00a3-4c82-93cd-bbe76493499a
Liatris cylindracea -91.28287 43.50806 idigbio 1899-07-29 a86b1039-3d65-4d56-9b8e-b21bae479005
Liatris cylindracea -93.16693 45.08003 idigbio NA ab08ecbd-0ae2-450a-a980-41bb722bd957
Liatris cylindracea -93.72596 45.46555 idigbio 1998-08-14 acfc182b-2c1a-4cd3-bf3e-95c3f195d957
Liatris cylindracea -90.88104 41.46936 idigbio 1892-08-28 adc2c1b0-2635-40f0-9e89-376767dbf0ac
Liatris cylindracea -82.75000 46.78333 idigbio 1981-08-19 af530cd6-8efd-4ec6-896b-dd1d869a73d6
Liatris cylindracea -92.77505 44.85476 idigbio 1985-06-02 afff1832-5d50-484b-995a-2f857a2851f8
Liatris cylindracea -92.38209 43.43849 idigbio 1978-07-28 b0202c66-fd93-4f67-83cd-3d4ea9e42cbf
Liatris cylindracea -83.72633 42.27087 idigbio 1861-07-31 b5b89b61-79b2-4c4b-84b9-d90b8771ac02
Liatris cylindracea -90.62689 41.55328 idigbio 1915-08-30 b9d3f69b-af0b-4dd4-afcc-4b89b6aa433d
Liatris cylindracea -83.72684 42.29218 idigbio NA ba55b9c5-0440-4675-b31e-8848c46a183e
Liatris cylindracea -87.59394 41.79420 idigbio 1893-08-30 c1e013cb-56e7-42df-bfea-6149d3c80be3
Liatris cylindracea -83.72633 42.27087 idigbio NA c43fd15a-fa5a-442e-855e-dba7e20ccbf0
Liatris cylindracea -87.59237 41.56148 idigbio 1963-08-15 c48ff884-3840-47d7-8e41-26f48e756986
Liatris cylindracea -92.81944 44.80028 idigbio 1988-07-26 c5d18b51-c0e0-41a5-89ff-18150389fdc4
Liatris cylindracea -92.01953 43.76751 idigbio 1994-07-29 cac2e14b-a141-46a2-8752-295077781dae
Liatris cylindracea -96.44879 47.73852 idigbio 1979-08-15 d2dffc26-23f6-4ada-80bd-089416868acb
Liatris cylindracea -92.74229 44.63663 idigbio 1993-08-04 d38599cb-2e1b-45ec-b36b-7de887b9dda9
Liatris cylindracea -83.74886 42.27913 idigbio NA d7d0aa52-634f-47d6-b075-943085e95051
Liatris cylindracea -91.36237 43.62304 idigbio 1966-07-29 d95fdaac-7492-4d9e-bae4-19697f5670fa
Liatris cylindracea -83.11906 42.68765 idigbio 1925-08-26 da2c61e9-70ae-4e10-a044-500e8c1da11c
Liatris cylindracea -91.63931 44.04996 idigbio 1888-09-13 dce4d4a7-953b-42a5-b090-2bb1c041b434
Liatris cylindracea -91.42269 43.75640 idigbio 1998-08-08 de5a1e0a-25cf-4516-b66d-7aae1d8f1422
Liatris cylindracea -93.28313 45.54365 idigbio NA e378dd93-1ad5-48b5-b6ae-ddd51ad69a21
Liatris cylindracea -93.13361 35.27833 idigbio 1972-10-14 e5f12f01-0eed-4e86-ae46-645186cd1978
Liatris cylindracea -93.28313 45.54365 idigbio NA e6ac07ba-e527-4b11-afc4-ef661015501f
Liatris cylindracea -92.66120 44.28943 idigbio NA e76aff04-082e-47a7-aa84-44ffa1195737
Liatris cylindracea -91.99041 44.17096 idigbio 1986-07-31 e9e8a887-f89a-4fdf-878e-9fd215d73c95
Liatris cylindracea -92.46120 38.14800 idigbio 2010-10-16 eabfcb0f-f8cc-4390-9b8c-c55208f45dd9
Liatris cylindracea -95.24502 45.44867 idigbio 1982-07-29 eb6a8b19-eaf3-4f3b-8184-00668eceb635
Liatris cylindracea -92.52110 44.56198 idigbio 1982-07-31 f7bc3b78-9217-4a1d-8901-a316455ac5ad
Liatris cylindracea -93.59000 44.02000 idigbio NA fc40ac35-8362-4ef4-be9c-434779d499f0
Liatris cylindracea -91.39316 43.95302 idigbio 1943-08-25 fca3412f-8420-43c8-b4b2-f6a18945d2b5
Liatris cylindracea -92.36995 44.41148 idigbio 1883-08-14 fdff27f2-68c0-4ec1-88f1-f1f3c177e009

Occurrence Data Cleaning

Modified based on script written by Charlotte Germain-Aubrey.
Modified and created by ML Gaynor.

Load Packages

library(dplyr)
library(tidyr)
library(rjson)  
library(RCurl) 
library(raster)
library(sp)

Load data files

The datafiles for this demo are stored in “data/cleaning_demo/”.

raw.data <- read.csv("data/cleaning_demo/SampleFL-data.csv")
head(raw.data)
##               species      lat      long year
## 1 Asclepias_curtissii 26.08294 -81.40069 1943
## 2 Asclepias_curtissii 28.29933 -70.70104 1955
## 3 Asclepias_curtissii 35.30096 -82.43857 1923
## 4 Asclepias_curtissii 26.14614 -80.45289  \\N
## 5 Asclepias_curtissii 26.57758 -81.92146 1972
## 6 Asclepias_curtissii 26.64141 -80.43702 1981

Exploring the datafile

How many species are included in this file?

(species_count <- raw.data %>%
                  group_by(species) %>%
                  tally())
## # A tibble: 5 x 2
##   species                   n
##   <fct>                 <int>
## 1 Asclepias_curtissii      38
## 2 Asimina_obovata          16
## 3 Pinus_australis          60
## 4 Pinus_palustris          14
## 5 Pityothamnus_obovatus    14

Taxonomic Cleaning

Why?
Lots of bad names due to misspelling, synonymns, and more.

How?
iPlant TNRS Tool.

In addition to accessing the iPlant TNRS Tool in R, you can also use the website to access the database. In this tutorial, we only focus on checking synoymns using one tool, however there are many other tools to check taxonomy. For example, taxize can be used to accesss taxomny databases including Encylopedia of Life (eol), Taxonomic Name Resolution Service (tnrs), Integrated Taxonomic Information Service (itis), and many more.

Using the iPlant TNRS tool

Transform the species names list

First, we need to subset the names from the list we made above.

# Subset from the list we made above
(names <- species_count$species)
## [1] Asclepias_curtissii   Asimina_obovata       Pinus_australis      
## [4] Pinus_palustris       Pityothamnus_obovatus
## 5 Levels: Asclepias_curtissii Asimina_obovata ... Pityothamnus_obovatus

Next, remove the underscores between the genus and species name, replace with spaces.

(names<-gsub('_',' ',names))
## [1] "Asclepias curtissii"   "Asimina obovata"       "Pinus australis"      
## [4] "Pinus palustris"       "Pityothamnus obovatus"

Turn the list into a string by adding commas.

(names<-paste(names,collapse=','))
## [1] "Asclepias curtissii,Asimina obovata,Pinus australis,Pinus palustris,Pityothamnus obovatus"

Using the package RCurl makes the string URL-encoded.

names<-curlEscape(names)
Query the API

Using the packages rjson and RCurl we can query the Application programming interface (API) for the iPlant TNRS tool.
We have to set an object equal to the API url address.

tnrs.api<-'http://tnrs.iplantc.org/tnrsm-svc'

Next, we send a request to the address specified above.

url<-paste(tnrs.api,'/matchNames?retrieve=best&names=',names,sep='')
tnrs.json<-getURL(url) 
Process the results

The response is saved as tnrs.json above. Now, the response needs to be converted to a readable format from JSON.

tnrs.results<-fromJSON(tnrs.json)

Next, we need to summarize the names and change the format to match our raw data. First we isolate the names submitted and the found accepted name. Next, we replace the space in between the genus and species names with an underscore. Then, we convert the data file to a dataframe. Finally, we rename the columns to match the names in raw.data.

corrected_names<-sapply(tnrs.results[[1]], function(x) c(x$nameSubmitted,x$acceptedName))
corrected_names <- gsub(" ", "_", corrected_names)
corrected_names<-as.data.frame(t(corrected_names),stringsAsFactors=FALSE)
names(corrected_names) <- c("species", "new")
corrected_names
##                 species                 new
## 1   Asclepias_curtissii Asclepias_curtissii
## 2       Asimina_obovata     Asimina_obovata
## 3       Pinus_australis     Pinus_palustris
## 4       Pinus_palustris     Pinus_palustris
## 5 Pityothamnus_obovatus     Asimina_obovata
Correcting names

We now need to correct the names in the raw data file. To do this, we are going to merge the datasets, so we can save the new name with all the information in our raw data file.

merged_datasets <- merge(raw.data, corrected_names, by = "species")
head(merged_datasets)
##               species      lat      long year                 new
## 1 Asclepias_curtissii 26.08294 -81.40069 1943 Asclepias_curtissii
## 2 Asclepias_curtissii 28.29933 -70.70104 1955 Asclepias_curtissii
## 3 Asclepias_curtissii 35.30096 -82.43857 1923 Asclepias_curtissii
## 4 Asclepias_curtissii 26.14614 -80.45289  \\N Asclepias_curtissii
## 5 Asclepias_curtissii 26.57758 -81.92146 1972 Asclepias_curtissii
## 6 Asclepias_curtissii 26.64141 -80.43702 1981 Asclepias_curtissii

Date Cleaning

In addition to taxonmy fixes, we need to check the format of other columns, like dates. If you noticed above, when the year is missing there is an odd symbol. Lets replace that symbol with “NA”

merged_datasets$year <- gsub("\\N", NA, merged_datasets$year)

Is there any other issues?
Visualize the year and data associated.

(year_count <- merged_datasets %>%
                  group_by(year) %>%
                  tally())
## # A tibble: 23 x 2
##    year      n
##    <chr> <int>
##  1 1896      1
##  2 1901      4
##  3 1923      6
##  4 1932      1
##  5 1943      1
##  6 1944      1
##  7 1955      3
##  8 1965      2
##  9 1972      2
## 10 1974      2
## # … with 13 more rows

Location Cleaning

Precision

How precise is the locality information?

Precision is influenced by the number of decimal places associated with latitude and longitude. See below how degree and distance are related.
places degree distance
0 1.0000 111 km
1 0.1000 11.1 km
2 0.0100 1.11 km
3 0.0010 111 m
4 0.0001 11.1 m

Here we are going to round to two decimal points.

(merged_datasets$lat <- round(merged_datasets$lat, digits = 2))
##   [1]  26.08  28.30  35.30  26.15  26.58  26.64  27.08  27.34  27.38  29.47
##  [11]  27.48  27.49  27.70  27.90  27.91  27.96  28.06  28.30  28.30  28.51
##  [21]  28.55  28.71  28.77  29.07  29.21  29.47  29.61  29.98  26.58  27.49
##  [31]  27.91  28.51  29.07  29.61   0.00  18.51 -28.55  28.71  28.71   0.00
##  [41]  28.85  29.07  28.71  29.21  29.28  29.47  29.61  29.68  29.98  26.95
##  [51]  27.91  28.51  29.28  29.98  28.30  28.51  28.55  28.71  28.71  28.00
##  [61]  28.77  29.07  29.21  29.28  29.47  29.59  29.61  29.68  29.73  29.80
##  [71]  29.90  29.91  29.95  29.98  29.99  30.02  30.04  30.15   0.00  30.23
##  [81]  30.23  30.24  30.32  30.33  30.41  30.42  30.45  30.46  30.49  30.50
##  [91]  30.58  30.61  30.61  30.61  30.61  30.67  30.70  30.80  30.87  29.61
## [101]  29.68  29.73  29.80  30.02  30.58  30.80  30.87  19.73  19.80  19.90
## [111]  19.91  18.95  19.98  19.99  26.58  26.90  26.95  27.19  27.34  27.39
## [121]  27.48  27.49  27.70  27.90  27.91  27.96  28.06  28.10  26.95  27.34
## [131]  27.38  27.48  27.70  27.91  27.96  28.06  28.30  28.30  28.51  28.55
## [141]  25.51  25.55
merged_datasets$long <- round(merged_datasets$long, digits = 2)

Removing impossible points

Prior to the next step we have 142 points, after we have 137.

merged_datasets <- merged_datasets %>%
                  filter(lat != 0, long != 0)

Remove points that are botanical gardens.

First, load the data file of the botanical gardens in florida.

bg.points <- read.csv("data/cleaning_demo/BotanicalGardensFloridaCoordinates.csv", header=TRUE)

Next, we are going to filter the merged_dataset to exclude any coordinates that are in the bg.points dataset.

merged_datasets <- merged_datasets %>%
            filter(!lat %in% bg.points$Lat & !long %in% bg.points$Long)

Removing duplicates

We currently have 137 observations. Once only unique records are obtained, we have 121 observation.

merged_datasets.unique <- merged_datasets %>%
                          distinct
head(merged_datasets.unique) 
##               species   lat   long year                 new
## 1 Asclepias_curtissii 26.08 -81.40 1943 Asclepias_curtissii
## 2 Asclepias_curtissii 28.30 -70.70 1955 Asclepias_curtissii
## 3 Asclepias_curtissii 35.30 -82.44 1923 Asclepias_curtissii
## 4 Asclepias_curtissii 26.58 -81.92 1972 Asclepias_curtissii
## 5 Asclepias_curtissii 26.64 -80.44 1981 Asclepias_curtissii
## 6 Asclepias_curtissii 27.08 -80.40 2000 Asclepias_curtissii

Trimming to the desired area

The desired area is Florida.

Using the package raster we are able to create a raster layer of the USA. If you want to download other countries using getData, check out the avaliable maps. Next, we can subset for Florida. State names are saved as “NAME_1” and county names are saved as “NAME_2” in the map file.

# name = 'GADM' which is Database of Global Administrative Areas
# level = 2 gives the  level of subdivision (States)
usa <- raster::getData('GADM', country='USA', level=2)

# subset
florida <- subset(usa, NAME_1=="Florida")

Next, convert merged dataset unique into a spatial file.

xy_data <- data.frame(x = merged_datasets.unique$long, y = merged_datasets.unique$lat)
coordinates(xy_data) <- ~ x + y
proj4string(xy_data) <- crs("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

Extract the values of the points for the Florida polygon. Then, recombine everything to gather all data. Next, filter for only the points in Florida, select the columns needed for MaxEnt input, and rename the new names as simply ‘species’.

over <- over(xy_data, florida)
total <- cbind(merged_datasets.unique, over)
florida_points <- total %>%
                  filter(NAME_1 == "Florida") %>%
                  dplyr::select(new, lat, long) %>%
                  rename(species = new)
                  #rename(replace = c("species" = "new")) # AEM EDIT
head(florida_points)
##               species   lat   long
## 1 Asclepias_curtissii 26.08 -81.40
## 4 Asclepias_curtissii 26.58 -81.92
## 5 Asclepias_curtissii 26.64 -80.44
## 6 Asclepias_curtissii 27.08 -80.40
## 7 Asclepias_curtissii 27.38 -80.44
## 8 Asclepias_curtissii 27.48 -82.36

Visualize the points using ggplot2.

library(ggplot2)
states <- map_data("state")
florida <- subset(states, region == "florida" )

ggplot()+ 
  geom_polygon(florida, mapping = aes(x = long, y = lat),color = "black", fill = "gray", size = .2) +
  geom_point(florida_points, mapping = aes(x = long, y = lat, col = species ))

Finally, write a csv for the maxent steps.

#We need a csv with only species, lat, long for MaxEnt input
write.csv(florida_points, "data/cleaning_demo/MaxEntPointsInput_cleaned.csv", row.names = FALSE)

Climate Processing

Original script by Charlotte Germain-Aubrey.
Modified and created by ML Gaynor.

Load Packages

library(maptools)
library(raster)
library(rgdal)
library(sp)
library(maps)
library(mapproj)

Load bioclim layers

alt_l <- raster("data/climate_processing/Bioclim/alt.bil")
bio1_l <- raster("data/climate_processing/Bioclim/bio1.bil")
bio2_l <- raster("data/climate_processing/Bioclim/bio2.bil")
bio3_l <- raster("data/climate_processing/Bioclim/bio3.bil")
bio4_l <- raster("data/climate_processing/Bioclim/bio4.bil")
bio5_l <- raster("data/climate_processing/Bioclim/bio5.bil")
bio6_l <- raster("data/climate_processing/Bioclim/bio6.bil")
bio7_l <- raster("data/climate_processing/Bioclim/bio7.bil")
bio8_l <- raster("data/climate_processing/Bioclim/bio8.bil")
bio9_l <- raster("data/climate_processing/Bioclim/bio9.bil")
bio10_l <- raster("data/climate_processing/Bioclim/bio10.bil")
bio11_l <- raster("data/climate_processing/Bioclim/bio11.bil")
bio12_l <- raster("data/climate_processing/Bioclim/bio12.bil")
bio13_l <- raster("data/climate_processing/Bioclim/bio13.bil")
bio14_l <- raster("data/climate_processing/Bioclim/bio14.bil")
bio15_l <- raster("data/climate_processing/Bioclim/bio15.bil")
bio16_l <- raster("data/climate_processing/Bioclim/bio16.bil")
bio17_l <- raster("data/climate_processing/Bioclim/bio17.bil")
bio18_l <- raster("data/climate_processing/Bioclim/bio18.bil")
bio19_l <- raster("data/climate_processing/Bioclim/bio19.bil")

Read map of desired extent

There are many ways to determine your desired extent -
for example, convex hulls are often reccomended. We reccomend you refer to ENM literature before selecting the extent for your own project.

FL <- rgdal::readOGR("data/climate_processing/FL/FLstate2.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/shellygaynor/Dropbox/Workshops/Botany2019WorkshopNotes/ENM Material - Gaynor/iDigBio_Rbased_Demo/data/climate_processing/FL/FLstate2.shp", layer: "FLstate2"
## with 200 features
## It has 13 fields

Crop bioclim layers to desired extent

Here is the starting layer.

plot(alt_l)

First, mask the bioclim layer with Florida. Notice that the plot is blank because the map is masked.

alt <- mask(alt_l, FL)
plot(alt)

Next, crop the bioclim layer to Florida’s extent.

alt <- crop(alt, extent(FL))
plot(alt)

Finally, write an .asc cropped file to be used for MaxEnt.

writeRaster(alt, "data/climate_processing/PresentLayers/alt.asc", format="ascii")

Repeat with all additional layers.

bio1 <- mask(bio1_l, FL)
bio1 <- crop(bio1, extent(FL))
#writeRaster(bio1, "data/climate_processing/PresentLayers/bio1.asc", format="ascii")

bio2 <- mask(bio2_l, FL)
bio2 <- crop(bio2, extent(FL))
#writeRaster(bio2, "data/climate_processing/PresentLayers/bio2.asc", format="ascii")

bio3 <- mask(bio3_l, FL)
bio3 <- crop(bio3, extent(FL))
#writeRaster(bio3, "data/climate_processing/PresentLayers/bio3.asc", format="ascii")

bio4 <- mask(bio4_l, FL)
bio4 <- crop(bio4, extent(FL))
#writeRaster(bio4, "data/climate_processing/PresentLayers/bio4.asc", format="ascii")

bio5 <- mask(bio5_l, FL)
bio5 <- crop(bio5, extent(FL))
#writeRaster(bio5, "data/climate_processing/PresentLayers/bio5.asc", format="ascii")

bio6 <- mask(bio6_l, FL)
bio6 <- crop(bio6, extent(FL))
#writeRaster(bio6, "data/climate_processing/PresentLayers/bio6.asc", format="ascii")

bio7 <- mask(bio7_l, FL)
bio7 <- crop(bio7, extent(FL))
#writeRaster(bio7, "data/climate_processing/PresentLayers/bio7.asc", format="ascii")

bio8 <- mask(bio8_l, FL)
bio8 <- crop(bio8, extent(FL))
#writeRaster(bio8, "data/climate_processing/PresentLayers/bio8.asc", format="ascii")

bio9 <- mask(bio9_l, FL)
bio9 <- crop(bio9, extent(FL))
#writeRaster(bio9, "data/climate_processing/PresentLayers/bio9.asc", format="ascii")

bio10 <- mask(bio10_l, FL)
bio10 <- crop(bio10, extent(FL))
#writeRaster(bio10, "data/climate_processing/PresentLayers/bio10.asc", format="ascii")

bio11 <- mask(bio11_l, FL)
bio11 <- crop(bio11, extent(FL))
#writeRaster(bio11, "data/climate_processing/PresentLayers/bio11.asc", format="ascii")

bio12 <- mask(bio12_l, FL)
bio12 <- crop(bio12, extent(FL))
#writeRaster(bio12, "data/climate_processing/PresentLayers/bio12.asc", format="ascii")

bio13 <- mask(bio13_l, FL)
bio13 <- crop(bio13, extent(FL))
#writeRaster(bio13, "data/climate_processing/PresentLayers/bio13.asc", format="ascii")

bio14 <- mask(bio14_l, FL)
bio14 <- crop(bio14, extent(FL))
#writeRaster(bio14, "data/climate_processing/PresentLayers/bio14.asc", format="ascii")

bio15 <- mask(bio15_l, FL)
bio15 <- crop(bio15, extent(FL))
#writeRaster(bio15, "data/climate_processing/PresentLayers/bio15.asc", format="ascii")

bio16 <- mask(bio16_l, FL)
bio16 <- crop(bio16, extent(FL))
#writeRaster(bio16, "data/climate_processing/PresentLayers/bio16.asc", format="ascii")

bio17 <- mask(bio17_l, FL)
bio17 <- crop(bio17, extent(FL))
#writeRaster(bio17, "data/climate_processing/PresentLayers/bio17.asc", format="ascii")

bio18 <- mask(bio18_l, FL)
bio18 <- crop(bio18, extent(FL))
#writeRaster(bio18, "data/climate_processing/PresentLayers/bio18.asc", format="ascii")

bio19 <- mask(bio19_l, FL)
bio19 <- crop(bio19, extent(FL))
#writeRaster(bio19, "data/climate_processing/PresentLayers/bio19.asc", format="ascii")

Selecting layers for MaxEnt

We only want to include layers that are not highly correlated. To assess which layers to include, we will look at the pearson correlation coefficient among layers.

#### Now let's do a correlation analysis to try and reduce the number of layers for MaxEnt
stack <- stack(bio1, bio2, bio3, bio4, bio5, bio6, bio7, bio8, bio9, bio10, bio11, bio12, bio13, bio14, bio15, bio16, bio17, bio18, bio19)
corr <- layerStats(stack, 'pearson', na.rm=TRUE)
c <- corr$`pearson correlation coefficient`

write.csv(c, "data/climate_processing/correlationBioclim.csv")

Option 1

Open this file in Excel to view the values. It is more or less up to you which layers are removed, as long as the resulting layer set is lowly correlated.

Option 2

If you do not want to decide which layers you should keep or remove, you can use the following tool. This tool is described as follows: “The absolute values of pair-wise correlations are considered. If two variables have a high correlation, the function looks at the mean absolute correlation of each variable and removes the variable with the largest mean absolute correlation.” - This shows you which layers to remove.

library(caret)
c <- data.matrix(read.csv("data/climate_processing/correlationBioclim.csv", header = TRUE, row.names = 1,
                           sep = ","))
c <- abs(c)
envtCor <- findCorrelation(c, cutoff = 0.80, names = TRUE, exact = TRUE)
sort(envtCor)
##  [1] "bio1"  "bio10" "bio11" "bio13" "bio15" "bio16" "bio17" "bio19" "bio4" 
## [10] "bio6"  "bio7"

Ecological Analysis using Points

This is based off scripts created by Hannah Owens and James Watling, and Anthony Melton.
Modified and created by ML Gaynor.

Load Packages

library(raster)
library(rgdal)
library(ENMTools)
library(dismo)
library(RStoolbox)
library(hypervolume)
library(ggplot2)
library(dplyr)
library(tidyr)
library(gtools)
library(ggbiplot)

Load Datafiles

Species occurrence records

This is the file you made at the end of the Occurrence_Data_Cleaning.R script.

florida_points <- read.csv("data/cleaning_demo/MaxEntPointsInput_cleaned.csv")

Subset occurrence records for each species

Asclepias_curtissii <- florida_points %>%
                       filter(species == "Asclepias_curtissii")

Asimina_obovata <- florida_points %>%
                   filter(species == "Asimina_obovata")

Pinus_palustris <- florida_points %>%
                   filter(species == "Pinus_palustris")

Raster layers

These files were made in the ClimateProcesssing.R script.

# Raster stack
list <- list.files("data/climate_processing/PresentLayers/", full.names = T, recursive = FALSE) 
list <- mixedsort(sort(list))
envtStack <- stack(list)

Removing highly correlated layers

This tool is described as follows: “The absolute values of pair-wise correlations are considered. If two variables have a high correlation, the function looks at the mean absolute correlation of each variable and removes the variable with the largest mean absolute correlation.”

The data matrix we use was made at the end of the ClimateProcesssing.R script.

library(caret)
c <- data.matrix(read.csv("data/climate_processing/correlationBioclim.csv", header = TRUE, row.names = 1,
                           sep = ","))
c <- abs(c)
envtCor <- findCorrelation(c, cutoff = 0.80, names = TRUE, exact = TRUE)
sort(envtCor)
##  [1] "bio1"  "bio10" "bio11" "bio13" "bio15" "bio16" "bio17" "bio19" "bio4" 
## [10] "bio6"  "bio7"

Subset uncorrelated layers

envt.subset<-subset(envtStack, c(3, 4, 6, 9, 10, 13, 15, 19)) 
envt.subset
## class      : RasterStack 
## dimensions : 155, 182, 28210, 8  (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : -87.625, -80.04167, 24.54167, 31  (xmin, xmax, ymin, ymax)
## crs        : NA 
## names      : bio2, bio3, bio5, bio8, bio9, bio12, bio14, bio18

ENM - Realized niche

For each occurrence record, extract the value for each bioclim variables using the package raster.

ptExtract_Asclepias_curtissii <- raster::extract(envt.subset, Asclepias_curtissii[3:2])
ptExtract_Asimina_obovata <- raster::extract(envt.subset, Asimina_obovata[3:2])
ptExtract_Pinus_palustris <- raster::extract(envt.subset, Pinus_palustris[3:2])

Convert the extracted information for all three species and combined. I wrote a quick function to convert the files to dataframes.

# Set-up the function
convert_ptExtract <- function(value, name){
                      # convert ptExtract to a data frame
                      value_df <- as.data.frame(value)
                      # add a column with species name
                      value_df_DONE <- value_df %>% 
                                      mutate(species = name)
                      # return data frame
                      return(value_df_DONE)
}

# Apply the function 
ptExtract_Asclepias_curtissii_df <- convert_ptExtract(ptExtract_Asclepias_curtissii, "Asclepias_curtissii")
ptExtract_Asimina_obovata_df <- convert_ptExtract(ptExtract_Asimina_obovata, "Asimina_obovata")
ptExtract_Pinus_palustris_df <- convert_ptExtract(ptExtract_Pinus_palustris, "Pinus_palustris")

# Combined the three converted files and remove any values where NA appears 
pointsamples_combined <- rbind(ptExtract_Asclepias_curtissii_df,ptExtract_Asimina_obovata_df, ptExtract_Pinus_palustris_df)
pointsamples_combined <- pointsamples_combined %>% 
                         drop_na(bio2, bio3, bio5, bio8, bio9, bio12, bio14, bio18)

PCA

Create two dataframes.

data.bioclim <- pointsamples_combined[, 1:8]
data.species <- pointsamples_combined[, 9]

Using only the bioclim columns to run the principal components analysis.

data.pca <- prcomp(data.bioclim, scale. = TRUE) 

Understanding the PCA - Optional

library(factoextra)
library(FactoMineR)

When you use the command prcomp your loading variables show up as rotational variables. Thanks to a really great answer on stack overflow you can even convert the rotational variable to show the relative contribution.

loadings <- data.pca$rotation
#summary(loadings)

There are two options to convert the loading to show the relative contribution, they both give the same answer so either can be used.

loadings_relative_A <- t(t(abs(loadings))/rowSums(t(abs(loadings))))*100
summary(loadings_relative_A)
##       PC1              PC2              PC3              PC4          
##  Min.   : 2.741   Min.   : 1.553   Min.   : 0.702   Min.   : 0.07019  
##  1st Qu.: 9.919   1st Qu.: 2.833   1st Qu.: 4.693   1st Qu.: 4.75169  
##  Median :13.826   Median : 5.577   Median : 8.274   Median :11.75173  
##  Mean   :12.500   Mean   :12.500   Mean   :12.500   Mean   :12.50000  
##  3rd Qu.:16.644   3rd Qu.:19.293   3rd Qu.:17.401   3rd Qu.:17.78491  
##  Max.   :17.173   Max.   :36.658   Max.   :31.438   Max.   :28.34724  
##       PC5               PC6              PC7              PC8        
##  Min.   : 0.2921   Min.   : 1.572   Min.   : 1.967   Min.   : 3.251  
##  1st Qu.: 5.1674   1st Qu.: 3.075   1st Qu.: 4.478   1st Qu.: 5.486  
##  Median :10.3572   Median : 6.785   Median : 8.171   Median : 8.745  
##  Mean   :12.5000   Mean   :12.500   Mean   :12.500   Mean   :12.500  
##  3rd Qu.:22.4759   3rd Qu.:18.127   3rd Qu.:20.815   3rd Qu.:17.589  
##  Max.   :25.4017   Max.   :36.447   Max.   :28.451   Max.   :29.207
loadings_relative_B <- sweep(x = abs(loadings), MARGIN = 2, STATS = colSums(abs(loadings)), FUN = "/")*100
summary(loadings_relative_B)
##       PC1              PC2              PC3              PC4          
##  Min.   : 2.741   Min.   : 1.553   Min.   : 0.702   Min.   : 0.07019  
##  1st Qu.: 9.919   1st Qu.: 2.833   1st Qu.: 4.693   1st Qu.: 4.75169  
##  Median :13.826   Median : 5.577   Median : 8.274   Median :11.75173  
##  Mean   :12.500   Mean   :12.500   Mean   :12.500   Mean   :12.50000  
##  3rd Qu.:16.644   3rd Qu.:19.293   3rd Qu.:17.401   3rd Qu.:17.78491  
##  Max.   :17.173   Max.   :36.658   Max.   :31.438   Max.   :28.34724  
##       PC5               PC6              PC7              PC8        
##  Min.   : 0.2921   Min.   : 1.572   Min.   : 1.967   Min.   : 3.251  
##  1st Qu.: 5.1674   1st Qu.: 3.075   1st Qu.: 4.478   1st Qu.: 5.486  
##  Median :10.3572   Median : 6.785   Median : 8.171   Median : 8.745  
##  Mean   :12.5000   Mean   :12.500   Mean   :12.500   Mean   :12.500  
##  3rd Qu.:22.4759   3rd Qu.:18.127   3rd Qu.:20.815   3rd Qu.:17.589  
##  Max.   :25.4017   Max.   :36.447   Max.   :28.451   Max.   :29.207

Plotting the PCA

First, I made a theme to change the background of the plot. Next, I changed the plot margins and the text size.

theme <- theme(panel.background = element_blank(),
               panel.border=element_rect(fill=NA),
               panel.grid.major = element_blank(),
               panel.grid.minor = element_blank(),
               strip.background=element_blank(),
               axis.ticks=element_line(colour="black"),
               plot.margin=unit(c(1,1,1,1),"line"), 
               axis.text = element_text(size = 12), 
               legend.text = element_text(size = 12), 
               legend.title = element_text(size = 12), 
               text = element_text(size = 12))

Next, ggbiplot where obs.scale indicates the scale factor to apply to observation, var.scale indicates the scale factor to apply to variables, ellipse as TRUE draws a normal data ellipse for each group, and circle as TRUE draws a correlation circle.

g <- ggbiplot(data.pca, obs.scale = 1, var.scale = 1, groups = data.species, ellipse = TRUE, circle = TRUE)
g <- g + theme
g <- g + scale_colour_manual(values = c("#73dfff", "#ff8700", "#7a8ef5"))
g <- g + theme(legend.direction = 'horizontal', legend.position = 'bottom')
g

***

Multivariate environment similarity surface (MESS) analysis and mahalanobis distances

For MESS, the more negative the score, the more likely you are to project into novel environments. For mahalanobis distances, zero equals completely overlapping. Comparatively, the higher the distance, the more likely you are to project into an environment that falls outside the range of the training region, leading to unreliable projections.

Read layer files

Training region layers
env.files.c <- list.files(path = "data/climate_processing/PresentLayers/", pattern = ".asc", full.names = TRUE) #Always check for tiff or asc
current <- stack(env.files.c)
names(current) <- c("alt", "bio1", "bio10", "bio11", "bio12", "bio13", "bio14", "bio15", "bio16", "bio17", "bio18", "bio19", "bio2", "bio3", "bio4", "bio5", "bio6", "bio7", "bio8", "bio9")
current <- setMinMax(current)
projection(current) <- "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs"
plot(current[[1]])

Projection layers

Projection layers can be a different geographic region or different time period.

env.files.p <- list.files(path = "data/climate_processing/Coastal_Plains_Layers/", pattern = ".asc", full.names = TRUE) #Always check for tiff or asc
proj <- stack(env.files.p)
names(proj) <- c("alt", "bio1", "bio10", "bio11", "bio12", "bio13", "bio14", "bio15", "bio16", "bio17", "bio18", "bio19", "bio2", "bio3", "bio4", "bio5", "bio6", "bio7", "bio8", "bio9")
proj <- setMinMax(proj)
projection(proj) <- "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs"
plot(proj[[1]])

MESS

For MESS, the more negative the score, the more likely you are to project into novel environments.
###### Generate data frame

ref <- as.data.frame(current) 
head(ref)
##   alt bio1 bio10 bio11 bio12 bio13 bio14 bio15 bio16 bio17 bio18 bio19 bio2
## 1  NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA   NA
## 2  87  188   264   104  1635   180    81    18   478   324   478   419  126
## 3  82  188   264   103  1627   179    81    18   476   323   476   417  126
## 4  87  187   264   103  1626   179    82    18   477   323   477   416  126
## 5  86  187   264   103  1625   179    82    18   478   323   478   415  127
## 6  78  187   263   103  1623   179    81    18   479   322   479   415  127
##   bio3 bio4 bio5 bio6 bio7 bio8 bio9
## 1   NA   NA   NA   NA   NA   NA   NA
## 2   42 6294  326   31  295  264  194
## 3   42 6304  326   31  295  264  194
## 4   42 6287  326   31  295  264  193
## 5   43 6291  325   31  294  264  193
## 6   43 6282  325   31  294  263  193
Run MESS
mess <- mess(x = proj, v = ref, full = FALSE) 
plot(mess)

Mahalanobis distance

For mahalanobis distances, zero equals completely overlapping. Comparatively, the higher the distance, the more likely you are to project into an environment that falls outside the range of the training region, leading to unreliable projections.

Greater numbers means the variables are more dissimilar.

Generate data frames

Generate data frames containing the environmental data.

refdat <- as.data.frame(ref) 
head(refdat)
##   alt bio1 bio10 bio11 bio12 bio13 bio14 bio15 bio16 bio17 bio18 bio19 bio2
## 1  NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA   NA
## 2  87  188   264   104  1635   180    81    18   478   324   478   419  126
## 3  82  188   264   103  1627   179    81    18   476   323   476   417  126
## 4  87  187   264   103  1626   179    82    18   477   323   477   416  126
## 5  86  187   264   103  1625   179    82    18   478   323   478   415  127
## 6  78  187   263   103  1623   179    81    18   479   322   479   415  127
##   bio3 bio4 bio5 bio6 bio7 bio8 bio9
## 1   NA   NA   NA   NA   NA   NA   NA
## 2   42 6294  326   31  295  264  194
## 3   42 6304  326   31  295  264  194
## 4   42 6287  326   31  295  264  193
## 5   43 6291  325   31  294  264  193
## 6   43 6282  325   31  294  263  193
projdat <- as.data.frame(proj) 
head(projdat)
##   alt bio1 bio10 bio11 bio12 bio13 bio14 bio15 bio16 bio17 bio18 bio19 bio2
## 1  NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA   NA
## 2  NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA   NA
## 3  NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA   NA
## 4  NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA   NA
## 5  NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA   NA
## 6  NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA   NA
##   bio3 bio4 bio5 bio6 bio7 bio8 bio9
## 1   NA   NA   NA   NA   NA   NA   NA
## 2   NA   NA   NA   NA   NA   NA   NA
## 3   NA   NA   NA   NA   NA   NA   NA
## 4   NA   NA   NA   NA   NA   NA   NA
## 5   NA   NA   NA   NA   NA   NA   NA
## 6   NA   NA   NA   NA   NA   NA   NA

Calculate the average and covariance matrix of the variables in the reference set.

ref.av <- colMeans(refdat, na.rm=TRUE)
ref.cov <- var(refdat, na.rm=TRUE)

Calculate the mahalanobis distance of each raster cell to the environmental center of the reference set for both the reference and the projection data set and calculate the ratio between the two.

mah.ref <- mahalanobis(x=refdat, center=ref.av, cov=ref.cov, tol=1e-20)
mah.pro <- mahalanobis(x=projdat, center=ref.av, cov=ref.cov, tol=1e-20)
mah.max <- max(mah.ref[is.finite(mah.ref)])
nt2 <- as.data.frame(mah.pro / mah.max)

Create and plot the raster layer.

NT2 <- read.asciigrid(env.files.p[[1]])
NT2@data <- nt2
NT2rast <- raster(NT2)
plot(NT2rast, col=rainbow(100))

***

Ecological Niche Modeling

MUST HAVE RJAVA INSTALLED
We will not demo these scripts during the workshop due to system requirements of rJava.
Script by Anthony Melton and ML Gaynor.

This script is for generating and testing ENMs using ENMEval. Please see https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12261 for the paper describing ENMEval and https://cran.r-project.org/web/packages/ENMeval/vignettes/ENMeval-vignette.html#cb2 for the vignette.

Set up java memory

options(java.parameters = "- Xmx16g") # increase memory that can be used

Load Packages

library(dplyr)
library(ENMeval)
library(rSDM)
library(dismo)
library(spThin)
library(readr)
library(rmaxent)
library(gtools)

Load Datafiles

Species occurrence records

This is the file you made at the end of the Occurrence_Data_Cleaning.R script.

florida_points <- read.csv("data/cleaning_demo/MaxEntPointsInput_cleaned.csv")

Subset occurrence records for each species

Asclepias_curtissii <- florida_points %>%
                       filter(species == "Asclepias_curtissii")

Asimina_obovata <- florida_points %>%
                   filter(species == "Asimina_obovata")

Pinus_palustris <- florida_points %>%
                   filter(species == "Pinus_palustris")

Raster layers

These files were made in the ClimateProcesssing.R script.

# Raster stack
list <- list.files("data/climate_processing/PresentLayers/", full.names = T, recursive = FALSE) 
list <- mixedsort(sort(list))
envtStack <- stack(list)

Removing highly correlated layers

This tool is described as follows: “The absolute values of pair-wise correlations are considered. If two variables have a high correlation, the function looks at the mean absolute correlation of each variable and removes the variable with the largest mean absolute correlation.”

The data matrix we use was made at the end of the ClimateProcesssing.R script.

library(caret)
c <- data.matrix(read.csv("data/climate_processing/correlationBioclim.csv", header = TRUE, row.names = 1,
                           sep = ","))
c <- abs(c)
envtCor <- findCorrelation(c, cutoff = 0.80, names = TRUE, exact = TRUE)
sort(envtCor)
##  [1] "bio1"  "bio10" "bio11" "bio13" "bio15" "bio16" "bio17" "bio19" "bio4" 
## [10] "bio6"  "bio7"

Subset uncorrelated layers

envt.subset<-subset(envtStack, c(3, 4, 6, 9, 10, 13, 15, 19)) 
envt.subset
## class      : RasterStack 
## dimensions : 155, 182, 28210, 8  (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : -87.625, -80.04167, 24.54167, 31  (xmin, xmax, ymin, ymax)
## crs        : NA 
## names      : bio2, bio3, bio5, bio8, bio9, bio12, bio14, bio18

Thin occurrence records

If there are still a lot of points after cleaning, you can further reduce to minimize biases in g-space.
The following code requires the thin_max.R script by Dan Warren available at http://enmtools.blogspot.com/.

nrow(Pinus_palustris)
## [1] 51
if (nrow(Pinus_palustris) > 50) {
source("functions/thin.max.R")
Pinus_palustris <- thin.max(Pinus_palustris, c("long", "lat"), 50) # trim down to N; check column names
}

nrow(Pinus_palustris)
## [1] 50
plot(x = Pinus_palustris$lat, y = Pinus_palustris$long, col = "black")

Designate background data

Randomly sample 10,000 background points from one background extent raster (only one per cell without replacement). Note: Since the raster has <10,000 pixels, you’ll get a warning and all pixels will be used for background. We will be sampling from the biome variable because it is missing some grid cells, and we are trying to avoid getting background points with NA.

bg <- randomPoints(envt.subset[[1]], n=10000)
bg <- as.data.frame(bg)

Check the bg point distribution in g-space

plot(envt.subset[[1]], legend=FALSE)

plot(x = bg$x, y = bg$y, pch = 16, col='red')

Model generation

ENMeval will generate multiple models and test them per the specified test method. Two important variables here are the regularization multiplier value and the feature class (fm). FC will allow for different shapes in response curves (linear, hinge, quadratic, product, and threshold) can be used in the model and RM will influence how many parameters are included in the model.

modeval <- ENMevaluate(occ = Pinus_palustris[, c("long", "lat")], 
                       env = envt.subset,
                       bg.coords = bg,
                       #categoricals = "Taxonomy",
                       algorithm = 'maxent.jar',
                       #RMvalues = c(0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4),
                       #fc = c("L", "H", "LQ", "LQH", "LQP", "LQPH", "LQPHT"),
                       RMvalues = c(1, 2.5, 5),
                       fc = c("L", "H", "LQH"), 
                       method = "block",
                       aggregation.factor = c(2,2),
                       #overlap = TRUE,
                       clamp = TRUE, 
                       rasterPreds = TRUE,
                       parallel = TRUE,
                       numCores = 4,
                       bin.output = TRUE,
                       progbar = TRUE)

Evaluating the Models

AICc score

Get the “best” model based on AICc score.

aic.opt <- modeval@models[[which(modeval@results$delta.AICc==0)]]
aic.opt
## class    : MaxEnt 
## variables: bio2 bio3 bio5 bio8 bio9 bio12 bio14 bio18 
## output html file no longer exists
aic.opt@lambdas
##  [1] "bio12, 0.0, 616.0, 1635.0"                    
##  [2] "bio14, 0.0, 23.0, 93.0"                       
##  [3] "bio18, 0.0, 236.0, 660.0"                     
##  [4] "bio2, 1.1480830825248387, 59.0, 138.0"        
##  [5] "bio3, -1.3121161746720424, 36.0, 54.0"        
##  [6] "bio5, 0.0, 310.0, 337.0"                      
##  [7] "bio8, 0.0, 257.0, 278.0"                      
##  [8] "bio9, 0.0, 123.0, 229.0"                      
##  [9] "linearPredictorNormalizer, 0.5513150681853487"
## [10] "densityNormalizer, 4864.835031469543"         
## [11] "numBackgroundPoints, 7807"                    
## [12] "entropy, 8.911469747354095"

Lambda

The Parse_lambdas.R script by John Baums can be found at https://github.com/johnbaums/rmaxent/blob/master/R/parse_lambdas.R. When features that have a lambda of > 0, these are the features that are included in the model.

source("functions/parse_lambdas.R")
parse_lambdas(aic.opt) 
## 
## Features with non-zero weights
## 
##  feature lambda min max   type
##     bio2  1.148  59 138 linear
##     bio3 -1.312  36  54 linear
## 
## 
## Features with zero weights
## 
##  feature lambda min  max   type
##    bio12      0 616 1635 linear
##    bio14      0  23   93 linear
##    bio18      0 236  660 linear
##     bio5      0 310  337 linear
##     bio8      0 257  278 linear
##     bio9      0 123  229 linear

Create a table with the stats for all models

results <- modeval@results
kable(results) %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", font_size = 10)) %>%
  scroll_box(width = "100%", height = "200px")
settings features rm train.AUC avg.test.AUC var.test.AUC avg.diff.AUC var.diff.AUC avg.test.orMTP var.test.orMTP avg.test.or10pct var.test.or10pct AICc delta.AICc w.AIC parameters AUC_bin.1 AUC_bin.2 AUC_bin.3 AUC_bin.4 diff.AUC_bin.1 diff.AUC_bin.2 diff.AUC_bin.3 diff.AUC_bin.4 test.or10pct_bin.1 test.or10pct_bin.2 test.or10pct_bin.3 test.or10pct_bin.4 test.orMTP_bin.1 test.orMTP_bin.2 test.orMTP_bin.3 test.orMTP_bin.4
L_1 L 1.0 0.6427 0.5950680 0.0769192 0.0951063 0.0145192 0.0833333 0.0277778 0.1233974 0.0395402 900.9929 9.8534940 0.0028874 8 0.5516992 0.4137206 0.8532727 0.5615794 0.1089778 0.1941263 0 0.0773210 0.0769231 0.4166667 0 0 0.0000000 0.3333333 0 0
H_1 H 1.0 0.7128 0.6340071 0.0881188 0.1274079 0.0330382 0.1233974 0.0395402 0.2019231 0.0605276 929.6835 38.5441091 0.0000000 18 0.5191248 0.4542398 0.9008484 0.6618153 0.2281117 0.2338202 0 0.0476998 0.3076923 0.5000000 0 0 0.0769231 0.4166667 0 0
LQH_1 LQH 1.0 0.7058 0.6245624 0.0879470 0.1276296 0.0286817 0.1041667 0.0434028 0.1442308 0.0575690 946.7588 55.6194488 0.0000000 21 0.5452257 0.4382712 0.9005281 0.6142244 0.1830345 0.2539500 0 0.0735339 0.0769231 0.5000000 0 0 0.0000000 0.4166667 0 0
L_2.5 L 2.5 0.6250 0.6162346 0.1108003 0.0983665 0.0297436 0.0625000 0.0156250 0.1410256 0.0282709 894.4267 3.2872847 0.0769727 4 0.4527741 0.4172217 0.8861131 0.7088297 0.2185759 0.1748903 0 0.0000000 0.2307692 0.3333333 0 0 0.0000000 0.2500000 0 0
H_2.5 H 2.5 0.6296 0.6047659 0.0595591 0.0894792 0.0256199 0.0416667 0.0069444 0.0817308 0.0138992 899.1838 8.0444577 0.0071339 6 0.4715196 0.4592780 0.7728074 0.7154583 0.2116216 0.1462953 0 0.0000000 0.0769231 0.2500000 0 0 0.0000000 0.1666667 0 0
LQH_2.5 LQH 2.5 0.6257 0.6278310 0.1056784 0.0860274 0.0243050 0.0817308 0.0138992 0.1410256 0.0282709 895.4963 4.3569047 0.0450891 5 0.4629622 0.4385594 0.8920545 0.7177480 0.2094973 0.1346121 0 0.0000000 0.2307692 0.3333333 0 0 0.0769231 0.2500000 0 0
L_5 L 5.0 0.6221 0.6222853 0.0983180 0.0806673 0.0293248 0.0400641 0.0021470 0.1394231 0.0264731 891.1394 0.0000000 0.3982554 2 0.4305357 0.4683350 0.8688504 0.7214199 0.2421766 0.0804926 0 0.0000000 0.3076923 0.2500000 0 0 0.0769231 0.0833333 0 0
H_5 H 5.0 0.6036 0.5674198 0.0217130 0.0529846 0.0252663 0.0000000 0.0000000 0.0000000 0.0000000 893.7779 2.6385224 0.1064667 2 0.4663615 0.5000000 0.6516588 0.6516588 0.2119382 0.0000000 0 0.0000000 0.0000000 0.0000000 0 0 0.0000000 0.0000000 0 0
LQH_5 LQH 5.0 0.6225 0.6171537 0.1022417 0.0855095 0.0342282 0.0400641 0.0021470 0.1586538 0.0365816 891.3237 0.1843087 0.3631947 2 0.4145146 0.4683350 0.8682198 0.7175452 0.2615454 0.0804926 0 0.0000000 0.3846154 0.2500000 0 0 0.0769231 0.0833333 0 0

Visualize model statistics

Relative occurrence rate
plot(modeval@predictions[[which(modeval@results$delta.AICc==0)]], main="Relative occurrence rate")

Variable importance
(df <- var.importance(aic.opt))
##   variable percent.contribution permutation.importance
## 1    bio12               0.0000                 0.0000
## 2    bio14               0.0000                 0.0000
## 3    bio18               0.0000                 0.0000
## 4     bio2              25.4036                45.1148
## 5     bio3              74.5964                54.8852
## 6     bio5               0.0000                 0.0000
## 7     bio8               0.0000                 0.0000
## 8     bio9               0.0000                 0.0000
barplot(df$permutation.importance, names.arg=df$variable, las=2, ylab="Permutation Importance")

delta.AICc
eval.plot(results)

Visualize Predicted Models

Plots all of the predictions by the different models

maps <- modeval@predictions
plot(maps)

Plots the predictions of the best model

plot(maps[[which(results$delta.AICc == 0)]], main = "Models with lowest AICc")

Response curves

Plots all of the response curves for the best model.

for (i in which(results$delta.AICc == 0)) {
  response(modeval@models[[i]])
}

The “best” model
# Pick "best" model based on AICc
which(results$delta.AICc == 0)
## [1] 7
n <- which(results$delta.AICc == 0)

# Isolate the "best" model
modeval@predictions[[n]]  # raw output; Check model name!
## class      : RasterLayer 
## dimensions : 155, 182, 28210  (nrow, ncol, ncell)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : -87.625, -80.04167, 24.54167, 31  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : L_5 
## values     : 0.00006318226, 0.0002055568  (min, max)
p <- predict(modeval@models[[n]], envt.subset) 
plot(p, zlim = c(0,1))  # plots the predictions in a 0-1 scale; logistic output

Save output

R saved datset

save(modeval, file = "PineENM.rda")
# load(file = "FILENAME.rda")

Raster

writeRaster(x = p, filename = "PineENM.asc", format = "ascii", NAFlag = "-9999", overwrite = T)



Alternative: ENMTools

Load Packages

library(raster)
library(ENMTools)
library(ENMeval)
library(gtools)
library(dplyr)
library(RStoolbox)
library(hypervolume)
library(phytools)
library(dismo)
#library(alphanull)

ENMTools

We are using ENMTools to generate ENM models using MAXENT. This package may cause you trouble, for troubleshooting, see set-up.

Correct data format to be used to generate models.

envt.subset <- setMinMax(envt.subset)
proj4string(envt.subset) <- crs("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

Correct species format to required ENMTools input.

Asclepias_curtissii_enm <- enmtools.species(species.name = "Asclepias_curtissii",
                                            presence.points = Asclepias_curtissii[,3:2])

Asclepias_curtissii_enm$range <- background.raster.buffer(Asclepias_curtissii_enm$presence.points, 25000, mask = envt.subset)
Asclepias_curtissii_enm$background.points = background.points.buffer(points = Asclepias_curtissii_enm$presence.points, radius = 5000, n = 10000, mask =  envt.subset[[1]])
##############################
Asimina_obovata_enm <- enmtools.species(species.name = "Asimina_obovata",
                                        presence.points = Asimina_obovata[,3:2])

Asimina_obovata_enm$range <- background.raster.buffer(Asimina_obovata_enm$presence.points, 25000, mask = envt.subset)
Asimina_obovata_enm$background.points = background.points.buffer(points = Asimina_obovata_enm$presence.points, radius = 5000, n = 10000, mask =  envt.subset[[1]])
#############################
Pinus_palustris_enm <- enmtools.species(species.name = "Pinus_palustris",
                                        presence.points = Pinus_palustris[,3:2])

Pinus_palustris_enm$range <- background.raster.buffer(Pinus_palustris_enm$presence.points, 25000, mask = envt.subset)
Pinus_palustris_enm$background.points = background.points.buffer(points = Pinus_palustris_enm$presence.points, radius = 5000, n = 10000, mask =  envt.subset[[1]])

Run each model

Use MAXENT and generate models. These models will be used in the next tab. Here, the proportion of data to withhold randomly is set to 0.25. To match the GUI settings, add rts.reps = 5.

Asclepias_curtissii_enm.mx <- enmtools.maxent(species = Asclepias_curtissii_enm, env = envt.subset, test.prop = 0.25)
Asimina_obovata_enm.mx <- enmtools.maxent(species = Asimina_obovata_enm, env = envt.subset, test.prop = 0.25)
Pinus_palustris_enm.mx <- enmtools.maxent(species = Pinus_palustris_enm, env = envt.subset, test.prop = 0.25)

Save models

save(Asclepias_curtissii_enm.mx, file = "data/enm_output/ENMTools/Asclepias_curtissii.rda")
save(Asimina_obovata_enm.mx, file = "data/enm_output/ENMTools/Asimina_obovata.rda" )
save(Pinus_palustris_enm.mx, file = "data/enm_output/ENMTools/Pinus_palustris.rda" )

Post-ENMs

Script by Anthony Melton and ML Gaynor.

Load Packages

library(raster)
library(ENMTools)
library(ENMeval)
library(gtools)
library(dplyr)
library(RStoolbox)
library(hypervolume)
library(phytools)
library(dismo)
#library(alphanull)

ENM models

You can also read the models you generated with the MaxEnt GUI into R.

Asclepias_curtissii_enm.mx.b <- raster("data/enm_output/Asclepias_curtissii/Asclepias_curtissii_avg.asc")
Asimina_obovata_enm.mx.b <- raster("data/enm_output/Asclepias_obovata/Asimina_obovata_avg.asc")
Pinus_palustris_enm.mx.b <- raster("data/enm_output/Pinus_palustris/Pinus_palustris_avg.asc")

See Ecological_Niche_Modeling.R to generate models in R.

load(file = "data/enm_output/ENMTools/Asclepias_curtissii.rda")
load(file = "data/enm_output/ENMTools/Asimina_obovata.rda")
load(file = "data/enm_output/ENMTools/Pinus_palustris.rda")

Bioclim layers

These steps are described in the previous script.

list <- list.files("data/climate_processing/PresentLayers/", full.names = T, recursive = FALSE) 
list <- mixedsort(sort(list))
envtStack <- stack(list)
envt.subset<-subset(envtStack, c(3, 4, 6, 9, 10, 13, 15, 19)) 
envt.subset <- setMinMax(envt.subset)
proj4string(envt.subset) <- crs("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

ENM Breadth

Niche breadth is the breadth of suitable climatic factors for a species, ranging from 0 to 1. When breadth is closer to 1 the more generalist species with wider climatic tolerances. Values closer to 0 represent more specialized species. The raster.breadth command in ENMTools measures the smoothness of suitability scores across a projected landscape. The higher the score, the more of the available niche space a species occupies.

Rbased

Ac_breadth <- ENMTools::raster.breadth(x = Asclepias_curtissii_enm.mx$suitability)
Ac_breadth$B2
## [1] 0.9884032
Ao_breadth <- ENMTools::raster.breadth(x = Asimina_obovata_enm.mx$suitability)
Ao_breadth$B2
## [1] 0.998853
Pp_breadth <- ENMTools::raster.breadth(x= Pinus_palustris_enm.mx$suitability)
Pp_breadth$B2
## [1] 0.9271576

MaxEnt GUI

Ac_breadth.b <- ENMTools::raster.breadth(x = Asclepias_curtissii_enm.mx.b)
Ac_breadth.b$B2
## [1] 0.7653468
Ao_breadth.b <- ENMTools::raster.breadth(x = Asimina_obovata_enm.mx.b)
Ao_breadth.b$B2
## [1] 0.9527736
Pp_breadth.b <- ENMTools::raster.breadth(x = Pinus_palustris_enm.mx.b)
Pp_breadth.b$B2
## [1] 0.8838691

ENM overlap

Calculating niche overlap, Schoener’s D, with ENMEval - Schoener’s D ranges from 0 to 1, where zero represents no similarity between the projections and one represents completely identical projections.

Rbased

### Stack the enm models and make sure the stack is named.
enm_stack <- stack(Asclepias_curtissii_enm.mx$suitability, Asimina_obovata_enm.mx$suitability,Pinus_palustris_enm.mx$suitability)
names(enm_stack) <- c("Asclepias_curtissii", "Asimina_obovata", "Pinus_palustris")

### Calculate 
calc.niche.overlap(enm_stack, stat = "D")
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |======================================================================| 100%
##                     Asclepias_curtissii Asimina_obovata Pinus_palustris
## Asclepias_curtissii                  NA              NA              NA
## Asimina_obovata               0.9677212              NA              NA
## Pinus_palustris               0.8795803       0.8852847              NA

MaxEnt GUI

## MaxEnt GUI files
### Stack the enm models and make sure the stack is named.
enm_stack.b <- stack(Asclepias_curtissii_enm.mx.b, Asimina_obovata_enm.mx.b,Pinus_palustris_enm.mx.b)
names(enm_stack.b) <- c("Asclepias_curtissii", "Asimina_obovata", "Pinus_palustris")

### Calculate 
calc.niche.overlap(enm_stack.b, stat = "D")
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |======================================================================| 100%
##                     Asclepias_curtissii Asimina_obovata Pinus_palustris
## Asclepias_curtissii                  NA              NA              NA
## Asimina_obovata               0.8156201              NA              NA
## Pinus_palustris               0.6427275       0.7844774              NA

Hypervolume

Niche space can be thought of as a multi-dimensional hypervolume. We’re using climatic data in this case, so we’re measuring the hypervolume of climatic niche space occupied by these species in FL.

The following command generates a binary predicted occurrence map. To run this command, a function must be generated. There are two functions in the Functions.R script. By sourcing the script, the functions will be generated and you can run them here!

source("functions/Functions_AEM.R") 

Load occurrence data

florida_points <- read.csv("data/cleaning_demo/MaxEntPointsInput_cleaned.csv")

Asclepias_curtissii.pts <- florida_points %>%
                           filter(species == "Asclepias_curtissii")

Asimina_obovata.pts <- florida_points %>%
                       filter(species == "Asimina_obovata")

Pinus_palustris.pts <- florida_points %>%
                       filter(species == "Pinus_palustris")

Create the binary maps

We will generate a binary map with any cell that contains a suitability score greater than or equal to the lowest score with an occurrence point as a presence. There are other ways to set the threshold, including percentiles or training statistics.

Asimina.dist <- make.binary.map(model = Asimina_obovata_enm.mx.b, occ.dat = Asimina_obovata.pts)
Pinus.dist <- make.binary.map(model = Pinus_palustris_enm.mx.b, occ.dat = Pinus_palustris.pts)
Asclepias.dist <- make.binary.map(model = Asclepias_curtissii_enm.mx.b, occ.dat = Asclepias_curtissii.pts)

Calculate hypervolume

Next, let’s work on getting some data from the predicted distributions! Niche space can be thought of as a multi-dimensional hypervolume. We’re using climatic data in this case, so we’re measuring the hypervolume of climatic niche space occupied by these species in FL.

Asimina.hv <- get_hypervolume(binary_projection = Asimina.dist, envt = envt.subset)
## 
## Building tree... 
## done.
## Ball query... 
## 
## done.
## Binding random points... done.
## Beginning volume calculation... done.
Pinus.hv <- get_hypervolume(binary_projection = Pinus.dist, envt = envt.subset)
## 
## Building tree... 
## done.
## Ball query... 
## 
## done.
## Binding random points... done.
## Beginning volume calculation... done.
Asclepias.hv <- get_hypervolume(binary_projection = Asclepias.dist, envt = envt.subset)
## 
## Building tree... 
## done.
## Ball query... 
## 
## done.
## Binding random points... done.
## Beginning volume calculation... done.

Compare the hypervolumes

hv_set <- hypervolume_set(hv1 = Asimina.hv, hv2 = Pinus.hv, check.memory = F)
## Choosing num.points.max=673639 (use a larger value for more accuracy.)
## Using minimum density of 1114.529403
## Retaining 111377 points in hv1 and 111584 points in hv2.
## Beginning ball queries... 
## 
## Building tree... 
## done.
## Ball query... 
## 
## done.
## 
## Building tree... 
## done.
## Ball query... 
## 
## done.
## Finished ball queries.
hypervolume_overlap_statistics(hv_set)
##       jaccard      sorensen frac_unique_1 frac_unique_2 
##     0.6414567     0.7815700     0.2177020     0.2191566
plot(hv_set)

get_volume(hv_set)
##                                              untitled 
##                                              99.93186 
##                                              untitled 
##                                             100.11801 
##                  Intersection of (untitled, untitled) 
##                                              78.17649 
##                         Union of (untitled, untitled) 
##                                             121.87338 
## Unique component of (untitled) relative to (untitled) 
##                                              21.75537 
## Unique component of (untitled) relative to (untitled) 
##                                              21.94152

Phylogenetic correlations

Let’s look at niche overlap over a tree!

Prepare data files

Load a tree file
tree <- read.newick(file = "data/AOC_Test_Demo/tree.tre")
tree <- ape::compute.brlen(phy = tree, method = "grafen")
plot(tree)

Add missing species
Liatris_chapmanii.pts <- read.csv("data/AOC_Test_Demo/clean_L_chapmaniiPointsSPOCC_EDITS.csv")[,2:3]
Generate ENM models
Asclepias_curtissii <- enmtools.species(species.name = "Asclepias_curtissii",
                                            presence.points = Asclepias_curtissii.pts[,3:2])
Asclepias_curtissii$range <- background.raster.buffer(Asclepias_curtissii$presence.points, 25000, mask = envt.subset)
Asclepias_curtissii$background.points = background.points.buffer(points = Asclepias_curtissii$presence.points, radius = 5000, n = 10000, mask =  envt.subset[[1]])
##############################
Asimina_obovata <- enmtools.species(species.name = "Asimina_obovata",
                                        presence.points = Asimina_obovata.pts[,3:2])
Asimina_obovata$range <- background.raster.buffer(Asimina_obovata$presence.points, 25000, mask = envt.subset)
Asimina_obovata$background.points = background.points.buffer(points = Asimina_obovata$presence.points, radius = 5000, n = 10000, mask =  envt.subset[[1]])
#############################
Pinus_palustris <- enmtools.species(species.name = "Pinus_palustris",
                                        presence.points = Pinus_palustris.pts[,3:2])
Pinus_palustris$range <- background.raster.buffer(Pinus_palustris$presence.points, 25000, mask = envt.subset)
Pinus_palustris$background.points = background.points.buffer(points = Pinus_palustris$presence.points, radius = 5000, n = 10000, mask =  envt.subset[[1]])
#############################
Liatris_chapmanii <- enmtools.species(species.name = "Liatris_chapmanii",
                                          presence.points = Liatris_chapmanii.pts)
Liatris_chapmanii$range <- background.raster.buffer(Liatris_chapmanii$presence.points, 25000, mask = envt.subset)
Liatris_chapmanii$background.points = background.points.buffer(points = Liatris_chapmanii$presence.points, radius = 5000, n = 10000, mask =  envt.subset[[1]])
Create “clade” object
clade=enmtools.clade(species = list(Asclepias_curtissii, Pinus_palustris, Asimina_obovata, Liatris_chapmanii), tree = tree)
check.clade(clade)
## 
## 
## An enmtools.clade object with 4 species
## 
## Species names: 
##   Asclepias_curtissii     Liatris_chapmanii   Asimina_obovata     Pinus_palustris
## 
## Tree: 
## 
## Phylogenetic tree with 4 tips and 3 internal nodes.
## 
## Tip labels:
## [1] "Asclepias_curtissii" "Liatris_chapmanii"   "Asimina_obovata"    
## [4] "Pinus_palustris"    
## 
## Rooted; includes branch lengths.
## 
## 
## Data Summary: 
## <table>
##  <thead>
##   <tr>
##    <th style="text-align:left;">   </th>
##    <th style="text-align:left;"> species.names </th>
##    <th style="text-align:left;"> in.tree </th>
##    <th style="text-align:left;"> presence </th>
##    <th style="text-align:left;"> background </th>
##    <th style="text-align:left;"> range </th>
##   </tr>
##  </thead>
## <tbody>
##   <tr>
##    <td style="text-align:left;"> Asclepias_curtissii </td>
##    <td style="text-align:left;"> Asclepias_curtissii </td>
##    <td style="text-align:left;"> TRUE </td>
##    <td style="text-align:left;"> 21 </td>
##    <td style="text-align:left;"> 92 </td>
##    <td style="text-align:left;"> present </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Liatris_chapmanii </td>
##    <td style="text-align:left;"> Liatris_chapmanii </td>
##    <td style="text-align:left;"> TRUE </td>
##    <td style="text-align:left;"> 84 </td>
##    <td style="text-align:left;"> 168 </td>
##    <td style="text-align:left;"> present </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Asimina_obovata </td>
##    <td style="text-align:left;"> Asimina_obovata </td>
##    <td style="text-align:left;"> TRUE </td>
##    <td style="text-align:left;"> 17 </td>
##    <td style="text-align:left;"> 69 </td>
##    <td style="text-align:left;"> present </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Pinus_palustris </td>
##    <td style="text-align:left;"> Pinus_palustris </td>
##    <td style="text-align:left;"> TRUE </td>
##    <td style="text-align:left;"> 51 </td>
##    <td style="text-align:left;"> 217 </td>
##    <td style="text-align:left;"> present </td>
##   </tr>
## </tbody>
## </table>

Correlation Tests

For range AOCs, an intercept >0.5 and negative slope are indicative of sympatric species, while an intercept of <0.5 and a positive slope are indicative on non-sympatric speciation. A low intercept and positive slope for niche overlap would indicate niche divergence.

Age-Range Correlation Test
range.aoc <- enmtools.aoc(clade = clade,  nreps = 10, overlap.source = "range")
summary(range.aoc)
## 
## 
## Age-Overlap Correlation test
## 
## 10 replicates 
## 
## p values:
##      (Intercept) empirical.df$age 
##        0.1818182        0.1818182

## NULL
Age-Overlap Correlation Test
glm.aoc <- enmtools.aoc(clade = clade,  env = envt.subset, nreps = 10, overlap.source = "glm")
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species Asclepias_curtissii 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species Liatris_chapmanii 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species Asimina_obovata 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species Pinus_palustris 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
summary(glm.aoc)
## 
## 
## Age-Overlap Correlation test
## 
## 10 replicates 
## 
## p values:
##      (Intercept) empirical.df$age 
##        0.7272727        0.9090909

## NULL

Phylogenetic diversity

Original script by J. Janzten. Modified by ML Gaynor.

Load Packages

library(raster)
library(ape)
library(phytools)
library(picante)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(rgdal)
library(rgeos)
library(sp)
library(spatialEco)

Load the tree file

tree <- read.tree(file = "data/phylogenetic_diversity/tree_withBL.tre")
tree <- ape::compute.brlen(phy = tree, method = "grafen")
plot(tree)

Read ENM models

Asclepias_curtissii_enm.mx.b <- raster("data/enm_output/Asclepias_curtissii/Asclepias_curtissii_avg.asc")
Asimina_obovata_enm.mx.b <- raster("data/enm_output/Asclepias_obovata/Asimina_obovata_avg.asc")
Pinus_palustris_enm.mx.b <- raster("data/enm_output/Pinus_palustris/Pinus_palustris_avg.asc")
Liatris_chapmanii_enm.mx.b <- raster("data/phylogenetic_diversity/Liatris_chapmanii/Liatris_chapmanii_avg.asc")

Reclassify and stack rasters

Reclassify

reclassify_raster <- function(OGraster){
  ## Reclassify the raster by the suitability score
  OGraster[OGraster >= 0.25] <- 1
  OGraster[OGraster < 0.25] <- 0
  OGraster[is.na(OGraster)] <- 0
  crs(OGraster) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
  return(OGraster)
}

Asclepias_curtissii_enm.mx.b  <- reclassify_raster(Asclepias_curtissii_enm.mx.b)
Asimina_obovata_enm.mx.b <- reclassify_raster(Asimina_obovata_enm.mx.b)
Pinus_palustris_enm.mx.b <- reclassify_raster(Pinus_palustris_enm.mx.b)
Liatris_chapmanii_enm.mx.b <- reclassify_raster(Liatris_chapmanii_enm.mx.b)

Inspect reclassified rasters

par(mfrow=c(1,4))
plot(Asclepias_curtissii_enm.mx.b)
plot(Asimina_obovata_enm.mx.b)
plot(Pinus_palustris_enm.mx.b)
plot(Liatris_chapmanii_enm.mx.b)

Stack and convert to dataframe

# Stack rasters and create dataframe
rasterstack <- stack(Asclepias_curtissii_enm.mx.b, Asimina_obovata_enm.mx.b,
                   Pinus_palustris_enm.mx.b, Liatris_chapmanii_enm.mx.b)
## Rename
names(rasterstack) <- c("Asclepias_curtissii1", "Asimina_obovata1", "Pinus_palustris1", "Liatris_chapmanii1")

## Convert to dataframe
rasterstack_df <- as.data.frame(rasterstack, xy = TRUE)

Extract ecoregions

Load Florida shapefile

Load the shapefile and correct the CRS (coordinate reference system), then crop to the area of interest.

## Load Florida shapefile
### Obtained at https://www.epa.gov/eco-research/ecoregion-download-files-state-region-4
shp <- ("data/phylogenetic_diversity/Florida_Raster/reg4_eco_l4/reg4_eco_l4.shp")
floridashape <- shapefile(shp)
## Correct CRS
newcrs <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
floridashape2 <- spTransform(floridashape, newcrs)
## Crop to florida
floridashape3 <- raster::crop(floridashape2, extent(-87.625, -80.04167, 24.54167, 31))
Convert dataframe to points
df.sp <- rasterstack_df[,1:2]
coordinates(df.sp)<-~x+y
crs(df.sp) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
Extract points from shapefile
df.sp.info <- sp::over(df.sp, floridashape3)
Bind dataframe with points and with ecosystems
rasterstack_df.info <- cbind(df.sp.info, rasterstack_df)
rasterstack_df.info2 <- rasterstack_df.info %>%
                        dplyr::group_by(L4_KEY) %>% 
                        dplyr::summarize(Asclepias_curtissii = max(Asclepias_curtissii1),
                               Asimina_obovata = max(Asimina_obovata1),
                               Pinus_palustris = max(Pinus_palustris1), 
                               Liatris_chapmanii = max(Liatris_chapmanii1))
rasterstack_df.info2df <- as.data.frame(rasterstack_df.info2, row.names = FALSE)
rasterstack_df.info4 <- rasterstack_df.info2df[1:20,]
rasterstack_df.info3 <- rasterstack_df.info4[,-1]
row.names(rasterstack_df.info3) <- NULL
row.names(rasterstack_df.info3) <- rasterstack_df.info4[,1]

kable(rasterstack_df.info3) %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", font_size = 10)) %>%
  scroll_box(width = "100%", height = "200px")
Asclepias_curtissii Asimina_obovata Pinus_palustris Liatris_chapmanii
65f Southern Pine Plains and Hills 0 1 1 1
65g Dougherty Plain 0 1 1 1
65h Tifton Upland 0 1 1 1
65o Tallahasee Hills/Valdosta Limesink 0 1 1 1
65p Southeastern Floodplains and Low Terraces 0 1 1 1
75a Gulf Coast Flatwoods 1 1 1 1
75b Southwestern Florida Flatwoods 1 1 1 1
75c Central Florida Ridges and Uplands 1 1 1 1
75d Eastern Florida Flatwoods 1 1 1 1
75e Okefenokee Plains 0 1 1 1
75f Sea Island Flatwoods 1 1 1 1
75g Okefenokee Swamp 0 1 1 1
75i Floodplains and Low Terraces 0 1 1 1
75j Sea Islands/Coastal Marsh 1 1 1 0
75k Gulf Barrier Islands and Coastal Marshes 1 1 1 1
75l Big Bend Coastal Marsh 1 1 1 1
76a Everglades 1 1 1 1
76b Big Cypress 1 1 1 1
76c Miami Ridge/Atlantic Coastal Strip 1 1 1 1
76d Southern Coast and Islands 1 1 1 1

Phylogenetic Diversity Calculations

Match tree and dataframe

matched <- match.phylo.comm(tree, rasterstack_df.info3)
# Make object for tree including matching taxa only (non-matching taxa will be pruned out)
matchtree <- matched$phy
# Make object for dataset including matching taxa only (non-matching taxa removed)
matchcomm <- matched$comm 

Calculate phylogenetic distance matrix for use in MPD and MNTD calculations

MPD is mean pairwise distance, or mean phylogenetic distance among all pairs of species within a community. MNTD is mean nearest taxon distance, or the mean distance between each species within a community and its closest relative.

phydist <- cophenetic(matched$phy)

Calculate indices

Null model options include taxa.labels (shuffle tips of phylogeny) amongothers Use ses.pd, ses.mpd and ses.mntd for more info on alternative null models and other parameters and output format. Number of runs - only using 99 runs (comparing to null model) due to time constraints. Ideally, you should have 1,000 - 5,000 runs.

PD

PD calculates Faith’s PD and standard effect size of PD. PD is phylogenetic distance.

pd_result <-ses.pd(matchcomm, matchtree, null.model="taxa.labels", runs=99)
MPD

Calculates mpd and ses mpd - equivalent to -NRI. NRI is net-relatedness index, standardizes MPD or mean pairwise distance. ses stands for standardized effect size.
If we had abundance data, we could include that for mpd and mntd. The default is abundance.weighted = FALSE.

mpd_result <- ses.mpd(matchcomm, phydist, null.model="taxa.labels", abundance.weighted = FALSE, runs=99)
MNTD

Calculates mntd and ses mntd - equivalent to -NTI. NTI or nearest taxon index, standardizes MNTD mean nearest taxon distance. ses stands for standardized effect size.

mntd_result <- ses.mntd(matchcomm, phydist, null.model="taxa.labels", runs=99)

Look at results

Positive ses values (obs.z values) and p values > 0.95 indicate overdispersion (greater than expected).
Negative ses values (obs.z values) and p values < 0.05 indicate clustering (less than expected).
Values not significantly different from zero indicate taxa in community randomly distributed across tree.

kable(pd_result) %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", font_size = 10)) %>%
  scroll_box(width = "100%", height = "200px")
ntaxa pd.obs pd.rand.mean pd.rand.sd pd.obs.rank pd.obs.z pd.obs.p runs
65f Southern Pine Plains and Hills 3 2.666667 2.447811 0.2658219 73.0 0.8233154 0.730 99
65g Dougherty Plain 3 2.666667 2.447811 0.2658219 73.0 0.8233154 0.730 99
65h Tifton Upland 3 2.666667 2.447811 0.2658219 73.0 0.8233154 0.730 99
65o Tallahasee Hills/Valdosta Limesink 3 2.666667 2.447811 0.2658219 73.0 0.8233154 0.730 99
65p Southeastern Floodplains and Low Terraces 3 2.666667 2.447811 0.2658219 73.0 0.8233154 0.730 99
75a Gulf Coast Flatwoods 4 3.000000 3.000000 0.0000000 50.5 NaN 0.505 99
75b Southwestern Florida Flatwoods 4 3.000000 3.000000 0.0000000 50.5 NaN 0.505 99
75c Central Florida Ridges and Uplands 4 3.000000 3.000000 0.0000000 50.5 NaN 0.505 99
75d Eastern Florida Flatwoods 4 3.000000 3.000000 0.0000000 50.5 NaN 0.505 99
75e Okefenokee Plains 3 2.666667 2.447811 0.2658219 73.0 0.8233154 0.730 99
75f Sea Island Flatwoods 4 3.000000 3.000000 0.0000000 50.5 NaN 0.505 99
75g Okefenokee Swamp 3 2.666667 2.447811 0.2658219 73.0 0.8233154 0.730 99
75i Floodplains and Low Terraces 3 2.666667 2.447811 0.2658219 73.0 0.8233154 0.730 99
75j Sea Islands/Coastal Marsh 3 2.666667 2.393939 0.2910752 76.0 0.9369648 0.760 99
75k Gulf Barrier Islands and Coastal Marshes 4 3.000000 3.000000 0.0000000 50.5 NaN 0.505 99
75l Big Bend Coastal Marsh 4 3.000000 3.000000 0.0000000 50.5 NaN 0.505 99
76a Everglades 4 3.000000 3.000000 0.0000000 50.5 NaN 0.505 99
76b Big Cypress 4 3.000000 3.000000 0.0000000 50.5 NaN 0.505 99
76c Miami Ridge/Atlantic Coastal Strip 4 3.000000 3.000000 0.0000000 50.5 NaN 0.505 99
76d Southern Coast and Islands 4 3.000000 3.000000 0.0000000 50.5 NaN 0.505 99
kable(mpd_result) %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", font_size = 10)) %>%
  scroll_box(width = "100%", height = "200px")
ntaxa mpd.obs mpd.rand.mean mpd.rand.sd mpd.obs.rank mpd.obs.z mpd.obs.p runs
65f Southern Pine Plains and Hills 3 1.777778 1.533109 0.2758344 78.0 0.8870138 0.780 99
65g Dougherty Plain 3 1.777778 1.533109 0.2758344 78.0 0.8870138 0.780 99
65h Tifton Upland 3 1.777778 1.533109 0.2758344 78.0 0.8870138 0.780 99
65o Tallahasee Hills/Valdosta Limesink 3 1.777778 1.533109 0.2758344 78.0 0.8870138 0.780 99
65p Southeastern Floodplains and Low Terraces 3 1.777778 1.533109 0.2758344 78.0 0.8870138 0.780 99
75a Gulf Coast Flatwoods 4 1.555556 1.555556 0.0000000 50.5 NaN 0.505 99
75b Southwestern Florida Flatwoods 4 1.555556 1.555556 0.0000000 50.5 NaN 0.505 99
75c Central Florida Ridges and Uplands 4 1.555556 1.555556 0.0000000 50.5 NaN 0.505 99
75d Eastern Florida Flatwoods 4 1.555556 1.555556 0.0000000 50.5 NaN 0.505 99
75e Okefenokee Plains 3 1.777778 1.533109 0.2758344 78.0 0.8870138 0.780 99
75f Sea Island Flatwoods 4 1.555556 1.555556 0.0000000 50.5 NaN 0.505 99
75g Okefenokee Swamp 3 1.777778 1.533109 0.2758344 78.0 0.8870138 0.780 99
75i Floodplains and Low Terraces 3 1.777778 1.533109 0.2758344 78.0 0.8870138 0.780 99
75j Sea Islands/Coastal Marsh 3 1.777778 1.616162 0.2350454 71.5 0.6875955 0.715 99
75k Gulf Barrier Islands and Coastal Marshes 4 1.555556 1.555556 0.0000000 50.5 NaN 0.505 99
75l Big Bend Coastal Marsh 4 1.555556 1.555556 0.0000000 50.5 NaN 0.505 99
76a Everglades 4 1.555556 1.555556 0.0000000 50.5 NaN 0.505 99
76b Big Cypress 4 1.555556 1.555556 0.0000000 50.5 NaN 0.505 99
76c Miami Ridge/Atlantic Coastal Strip 4 1.555556 1.555556 0.0000000 50.5 NaN 0.505 99
76d Southern Coast and Islands 4 1.555556 1.555556 0.0000000 50.5 NaN 0.505 99
kable(mntd_result) %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", font_size = 10)) %>%
  scroll_box(width = "100%", height = "200px")
ntaxa mntd.obs mntd.rand.mean mntd.rand.sd mntd.obs.rank mntd.obs.z mntd.obs.p runs
65f Southern Pine Plains and Hills 3 1.555556 1.299663 0.2863549 73.5 0.8936193 0.735 99
65g Dougherty Plain 3 1.555556 1.299663 0.2863549 73.5 0.8936193 0.735 99
65h Tifton Upland 3 1.555556 1.299663 0.2863549 73.5 0.8936193 0.735 99
65o Tallahasee Hills/Valdosta Limesink 3 1.555556 1.299663 0.2863549 73.5 0.8936193 0.735 99
65p Southeastern Floodplains and Low Terraces 3 1.555556 1.299663 0.2863549 73.5 0.8936193 0.735 99
75a Gulf Coast Flatwoods 4 1.166667 1.166667 0.0000000 50.5 NaN 0.505 99
75b Southwestern Florida Flatwoods 4 1.166667 1.166667 0.0000000 50.5 NaN 0.505 99
75c Central Florida Ridges and Uplands 4 1.166667 1.166667 0.0000000 50.5 NaN 0.505 99
75d Eastern Florida Flatwoods 4 1.166667 1.166667 0.0000000 50.5 NaN 0.505 99
75e Okefenokee Plains 3 1.555556 1.299663 0.2863549 73.5 0.8936193 0.735 99
75f Sea Island Flatwoods 4 1.166667 1.166667 0.0000000 50.5 NaN 0.505 99
75g Okefenokee Swamp 3 1.555556 1.299663 0.2863549 73.5 0.8936193 0.735 99
75i Floodplains and Low Terraces 3 1.555556 1.299663 0.2863549 73.5 0.8936193 0.735 99
75j Sea Islands/Coastal Marsh 3 1.555556 1.292929 0.2881093 74.0 0.9115510 0.740 99
75k Gulf Barrier Islands and Coastal Marshes 4 1.166667 1.166667 0.0000000 50.5 NaN 0.505 99
75l Big Bend Coastal Marsh 4 1.166667 1.166667 0.0000000 50.5 NaN 0.505 99
76a Everglades 4 1.166667 1.166667 0.0000000 50.5 NaN 0.505 99
76b Big Cypress 4 1.166667 1.166667 0.0000000 50.5 NaN 0.505 99
76c Miami Ridge/Atlantic Coastal Strip 4 1.166667 1.166667 0.0000000 50.5 NaN 0.505 99
76d Southern Coast and Islands 4 1.166667 1.166667 0.0000000 50.5 NaN 0.505 99
Plot PD results by community
ggplot()+
  geom_point(data = pd_result, aes(x = rownames(pd_result), y = pd.obs))+
  ggtitle("PD values for FL by ecoregion")+
  xlab("Community")+
  ylab("PD")+
  theme(text = element_text(size=10), axis.text.x = element_text(angle = 90, hjust = 1))

Plot MPD results by community
ggplot()+
  geom_point(data = mpd_result, aes(x = rownames(mpd_result), y = mpd.obs))+
  ggtitle("MPD values for FL by ecoregion")+
  xlab("Community")+
  ylab("MPD")+
  theme(text = element_text(size=10), axis.text.x = element_text(angle = 90, hjust = 1))

Plot MNTD results by community
ggplot()+
  geom_point(data = mntd_result, aes(x = rownames(mntd_result), y = mntd.obs))+
  ggtitle("MNTD values for FL by ecoregion")+
  xlab("Community")+
  ylab("MNTD")+
  theme(text = element_text(size=10), axis.text.x = element_text(angle = 90, hjust = 1))